# Human Motion Estimation Based on Low Dimensional Space Incremental Learning

Authors: Wanyi Li; Jifeng Sun
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

This paper proposes a novel algorithm called low dimensional space incremental learning (LDSIL) to estimate the human motion in 3D from the silhouettes of human motion multiview images. The proposed algorithm takes the advantage of stochastic extremum memory adaptive searching (SEMAS) and incremental probabilistic dimension reduction model (IPDRM) to collect new high dimensional data samples. The high dimensional data samples can be selected to update the mapping from low dimensional space to high dimensional space, so that incremental learning can be achieved to estimate human motion from small amount of samples. Compared with three traditional algorithms, the proposed algorithm can make human motion estimation achieve a good performance in disambiguating silhouettes, overcoming the transient occlusion, and reducing estimation error.


---

## Body

## 1. Introduction

Human motion estimation has become a hot researching topic [1–3], but it is a challenging task. Unusually, we are very interested in estimating the human motion in 3D from the silhouettes of human motion multiview images. The challenges are mainly as follows: firstly, it is hard to build the mapping between multiview silhouettes and human motion in 3D; secondly, the matching is ambiguous between multiview silhouettes and human motion in 3D; finally, it is hard to determine the spatial position information of the human motion depicted in the multiview images. In the past few years, a number of algorithms have been proposed to estimate the human motion. In the works of Sigal et al. [4] and Deutscher and Reid [5], they use some improved particle filters to estimate the human motion in 3D. It cannot work well because of searching (sampling) in high dimensional (HD) space for many times. There are some problems if we directly search in HD space, for example, searching in large scale will get the invalid data and searching in small scale will not get the target data. Moreover, if searching many times in small scale, it also brings about the invalid data. Li et al. use principal component analysis (PCA) to reduce the dimension of the HD data samples and simulated annealing particle swarm optimism (SAPSO) to estimate the human motion in 3D [6]. This algorithm is time-consuming and its performance is not so well because the HD data converted from corresponding low dimensional (LD) data will be quite different from original HD data. Besides, it does not consider the spatial position information of the human motion. Some traditional Monte Carlo methods [7–9] have drawbacks, which can not ensure collecting the best sample each time during searching, thus stochastic extremum memory adaptive searching (SEMAS) is proposed to solve this problem. In the work of Wang et al. [10], Gaussian process dynamical model (GPDM) can be used to reduce the dimension of the HD data to acquire the corresponding LD data and build the mapping from LD space to HD space, but GPDM can not quickly reduce the dimension of the new HD data sample and acquire the new corresponding LD data; thus incremental probabilistic dimension reduction model (IPDRM) is proposed to solve this problem based on GPDM. Some improved incremental or nonincremental learning algorithms in [11–14] can not satisfy our need. The limitations are that output data denoting the class label or other simple information have only one dimension, which can not satisfy the description of some output data, and they can not carry out the unsupervised incremental learning of HD data. Inspired by the researches stated above, the key to estimate the human motion in 3D depends on generating the better prior information. The human motion in 3D can be estimated more accurately through searching around the better prior information in small scale only once. In this paper, we mainly focus on the regular human motion cycle (walking or running).Our task is how to use the small amount of samples of HD data as the prior information to estimate the human motion in 3D which matches the multiview images. Based on the works of the researches above, we proposed a novel algorithm called low dimensional space incremental learning (LDSIL). The LDSIL mainly carries out through SEMAS and IPDRM to collect the new HD samples and updates the mapping from the LD space to HD space through the selection of new HD samples, thus the searching in the LD space can generate the more accurate HD data to estimate the corresponding human motion in 3D. Then, SEMAS is used to find the spatial position of human motion in 3D, and it can find the best data sample during searching more easily. IPDRM is used to reduce the dimension of the new HD data sample and acquire the new corresponding LD data, and it can help to select the new HD samples through the mapping of incremental dimension reduction. Moreover, it provides the LD space to generate the valid HD data. Based on IPDRM, the method of selecting the HD data samples for incremental learning can be achieved by comparing corresponding LD data.The main contribution of this paper is listed as follows:(1)
SEMAS is proposed to find the spatial position of human motion model. It can get the best data sample better than the traditional Monte Carlo methods.(2)
IPDRM is proposed to reduce the dimension of the new HD data sample and acquire the new corresponding LD data. It can promote the incremental learning in LD space to update the mapping from the LD space to HD space. Besides, it provides the LD space to generate valid HD data through searching.(3)
The method of selecting the HD data samples is proposed; it is used to update the mapping from LD space to HD space.Overall, due to LDSIL being able to make use of LD data, it can solve these problems mentioned above and contribute a lot to estimating the human motion in 3D, which has the better performance than other traditional algorithms, including disambiguating silhouettes, overcoming the transient occlusion, and reducing estimation error.The rest of this paper is organized as follows: Section2 introduces the corresponding data and models, and they are used to estimate the human motion. Section 3 proposes SEMAS algorithm to find the spatial position of human motion model. Section 4 proposes IPDRM to achieve the incremental dimension reduction. Section 5 proposes orthogonal least squares learning of multiple outputs (OLSLMO) to learn the mapping (HD space to LD space and LD space to HD space). Section 6 discusses LDSIL based on Sections 2–5: this section mentions how to select the new HD data sample from the estimated human motion models to achieve the incremental learning in the LD space. Section 7 proposes the method of searching in LD space to estimate the human motion model; the method is taking the advantage of the SEMAS and IPDRM. Section 8 shows the validity of proposed algorithm (LDSIL) through the experiments and evaluations. Section 9 discusses the limitation of LDSIL algorithm and the improvement in the future.Let us give the more detailed discussion in the following sections.


## 2. Corresponding Data and Models

We introduce the corresponding data and models in the works [4, 5]. All image data can be found in HumanEva-I dataset [4], as shown in Figures 1(a)–1(d). Figure 1(a) shows the human motion model denoting the human motion in 3D, which is described by HD data. The model is our estimated object, which need match the limbs in the multiview images. Figure 1(b) shows the multiview images, which are depicting human motion and its spatial position. After using some image segmentation algorithms [15–17] to process the multiview images, we can get the silhouettes as shown in Figure 1(c). Then, we project the model to the corresponding views and obtain the projection images in Figure 1(d). The images in Figure 1(d) are used to compare with the images in Figure 1(c).Figure 1
Data and models.
(a)
Human motion model
(b)
Multiview images
(c)
Silhouette images
(d)
Projection images
(e)
Corresponding LD data of HD data samples (3D)In the following, we give some definitions for the abovementioned data and models. The equation is built as below:(1)wti=pmt∣ypti=∏k=1Nvpmtk∣ypti=exp⁡-SNv∑k=1Nv1-αBtkBtk+Htk+αRtkRtk+Htk,where t denotes the tth frame, k denotes the kth view of camera, Nv is the number of views, α denotes the weight coefficient, p denotes the pixel of the image, Mtf denotes the set of pixels in the silhouette image, and Mtb denotes the set of pixels in the projection image of the human motion model. Then, Mtf(p) is statistical function of p, which is in the silhouette (Mtf(p)=1, p∈Mtf or Mtf(p)=0, p∉Mtf). Mtb(p) is statistical function of p, which is in the model projection (Mtb(p)=1, p∈Mtb or Mtb(p)=0, p∉Mtb). Thus, in the kth view, Rtk=∑p(Mtf(p)(1-Mtb(p))) denotes the pixel number of p in the silhouette and not in the model projection, Btk=∑p(Mtb(p)(1-Mtf(p))) denotes the pixel number of p in the model projection and not in the silhouette, and Htk=∑p(Mtf(p)Mtb(p)) denotes the pixel number of p in the silhouette and the model projection. Then, the object image consists of multiview silhouette images, mtk denotes the kth view feature of object image, mt denotes multiview features of object image, ypti∈RD1 and yupti∈RD2, respectively, denote the HD data samples which contain spatial position information (6 dimensions in total) and no spatial position information, D1-D2=6, and thus p(mtk∣ypti) is the conditional probability of matching image feature of the kth view from the appearance of ypti. Usually, we let S=100 and α=0.5 in most cases. After dimension reduction, the corresponding LD data of the HD data samples with no spatial position information can be seen in Figure 1(e). The HD data samples are the small amount of samples which cover one cycle of human motion.For the better description of our proposed algorithm in the following section, we need to define the symbols of some operations as follows.yq1:q2 is extracting the q1th-q2th dimension data of vector y as a subvector; A1 and B1 are matrices or vectors, A1*B1=D1, Di,j1=A1i,jB1i,j; ones(n1,n2)∈Rn1×n2 is the matrix whose elements are 1; rand(n1,n2)∈Rn1×n2 is the matrix whose elements are conforming the distribution U(-1,1); the set can be described by {m1,m1+1,…,m2}=m1:m2, m1<m2, m1,m2∈N; H=H2=(∑i∑jHi,j2)1/2 is the calculation of 2-norm (Euclidean-norm), and H can denote the vector or matrix.


## 3. Stochastic Extremum Memory Adaptive Searching

We will propose SEMAS algorithm in this section, which is used to find the spatial position of human motion model. The spatial position information includes the angles and coordinates of root marker in the model, which is the low dimensional data. In timet (the tth frame), let yptiq1:q2 denote the spatial position information of ypti∈RD1, the set Cq={q1,q1+1,…,q2}=q1:q2, Cq∪C-q={1,2,…,D1}, and ypi=yptiq1:q2 and ya=ypti(C-q); we get(2)wti=pmt∣ypti=pmt∣yptiq1:q2,yptiC-q=pmt∣ypi,ya=fpypi,where the SEMAS can be carried out as follows. We denote k⌢1 as the ordinal number of maximum element in Wo={wt1,wt2,…,wtNe}, wtk⌢1>wti, i=1,2,…,Ne, i≠k⌢1, and yrpi is a sample of the rth time searching in time t, yrpi∈Rd, d=q2-q1+1, 1≤d≤6. Then, we let uti=v*rand(d,1), v∈Rd is the scale parameter vector, ηi=rand(1,1), and L is the searching parameter of the step size; thus the Ne-1 searching vectors can be given by u~ti=ηiL(uti/uti) or Luti, i=1,2,…,Ne-1; the Ne-1 new data samples can be obtained by yr+1pi=yrpk⌢1+u~ti, i=1,2,…,Ne-1. Here, the extremum needs to be kept, let yr+1pNe=yrpk⌢1, and, in the next time searching, it is recording the extremum for ensuring obtaining the data samples which is not worse than the ones in the last time searching. In each searching, we can get Ne data samples. After searching twice, we can adjust the value of L according to the equation as below:(3)L=A1+γ1yr-1pk⌢1-yr-2pk⌢1,pmt∣yr-1pk⌢1,ya-pmt∣yr-2pk⌢1,ya<η′,A2+γ2yr-1pk⌢1-yr-2pk⌢1,pmt∣yr-1pk⌢1,ya-pmt∣yr-2pk⌢1,ya≥η′,r>1,γ1≥γ2≥0,A1≥A2>0,η′>0,where when the subtraction of the best object values in the past two times searching is below some value, L can be enlarged to search better data samples; on the contrary, if the subtraction is above some value, it illustrates that the value of L can be shrunk to avoid missing the best data sample. Thus, the method can make the convergence faster. After r times searching (r is large enough), when fp(ypi) is unimodal approximately, the best sample yp* can be obtained by(4)yp*=argmax⁡ypifpypi≈yrpk⌢1,and the derivation of (4) and the pseudocode can be seen in the appendices.Then, we compare and analyze the performances between traditional Monte Carlo method and SEMAS as depicted in Figures2(a) and 2(b). In Figure 2, we take 4 times searching and 3 data samples each time for example. Figure 2(a) shows that the traditional Monte Carlo method can not ensure that the searched data samples are not worse than the last time, which keeps the larger weight data samples simple. The reason is that it does not keep the best data sample each time during searching, which can not be compared with the data samples searched in the next time; moreover, the length and direction of searching vector are stochastic and uncontrolled. After several times searching, in most cases, it may appear that all of the searched data samples can not close the fp*=max⁡(fp) better, and the mean f-p=∑w~tiyrpi,  w~ti=fp(yrpi)/∑fp(yrpi) can not reach fp* also. As shown in Figure 2(b), the SEMAS will keep the best data sample in each time searching, which will be compared with the searched data samples in the next time searching. It can ensure that the searched data samples can not be worse than the last time searching and can adjust the length of searching vector to control the searching scale, according to subtraction of the two best values of objective function fp in the past twice searching; thus it will close the fp*=max⁡(fp) better and has more chance to get the better data samples than traditional Monte Carlo method.Figure 2
The comparison of traditional Monte Carlo method and SEMAS.
(a)
Traditional Monte Carlo method
(b)
SEMAS


## 4. Incremental Probabilistic Dimension Reduction Model

### 4.1. Probabilistic Models of GPDM

We give the probabilistic models of the GPDM in [10] as follows:(5)pY∣X,β¯,W=WN2πNDKYDexp⁡-12tr⁡KY-1YW2YT,(6)pX∣α-=px12πN-1qKxqexp⁡-12tr⁡KX-1X2:NX2:NT,In (5), HD data sequence can be denoted by Y=y1,…,yi,…,yNT∈RN×D, yi∈RD, and LD data sequence can be denoted by X=[x1,…,xi,…,xN]T∈RN×q, xi∈Rq. Kernel matrix with parameter β¯=[β1,β2] is KY∈RN×N, KYi,j=kYxi,xj=exp⁡(-β1/2xi-xj2)+β2-1δxixj; W is scale parameter matrix, p(W)=∏m=1D(2/κ2π)exp⁡(-wm2/2κ2), wm>0, κ=10-3; in (6), kernel matrix with parameter α-=[α1,α2,α3,α4] is KX∈RN-1×N-1, (KX)i,j=kX(xi,xj)=α1exp⁡(-α2/2xi-xj2)+α3xiTxj+α4-1δxixj, X1:N-1=[x1,x2,…,xN-1]T, X2:N=[x2,x3,…,xN]T, and x1 conforms Gaussian distribution of q dimensions. α- and β¯ satisfy p(α¯)∝∏iαi-1 and p(β¯)∝∏iβi-1, respectively. Then, we give the conditional distribution as follows:(7)pY*∣X*,Γ=WM2πMDKY*Dexp⁡-12tr⁡KY*-1ZYW2ZYT,(8)pX*∣Γ=px1*2πM-1qKX*qexp⁡-12tr⁡KX*-1ZXZXT,where Γ={X,Y,α-,β-,W}, Y*∈RM×D, X*∈RM×q, ZY=Y*-ATKY-1Y, KY*=B-ATKY-1A, ZX=X2:N*-CTKX-1X2:N, and KX*=D-CTKX-1C. Kernel matrices include A∈RN×M, Ai,j=kY(xi,xj*), B∈RM×M, Bi,j=kY(xi*,xj*), C∈R(N-1)×(M-1), Ci,j=kX(xi,xj*) and D∈R(M-1)×(M-1), and Di,j=kX(xi*,xj*).When training GPDM,Y is known; thus p(Y) is constant, and min⁡(-ln⁡p(X,α¯,β¯,W∣Y))⇔min⁡(-ln⁡p(X,α¯,β¯,W,Y)). Then, the LD data and corresponding parameters can be calculated as below:(9)X,α¯,β¯,W=argmin⁡X,α¯,β¯,W(-ln⁡p(X,α¯,β¯,W,Y))=argmin⁡X,α¯,β¯,W-ln⁡pY∣X,β¯,WpX∣α-pα¯pβ¯pW·pX∣α-pα¯pβ¯pW.Optimization of (9) can use the method of scaled conjugate gradient (SCG) [18, 19]. After the optimization, X can be depicted as shown in Figure 1(e). According to (7), y*~N(μY(x*),σY2(x*)I), the mapping from LD data to HD data can be given by the mean as below:(10)y*=fx*=μY(x*)=YTKY-1kYx1,x*,kYx2,x*,…,kYxN,x*T=YTKY-1kYx*.


### 4.2. The Mapping of Incremental Dimension Reduction

After finishing training, GPDM can not process the new HD data sample, which means it can not acquire the corresponding LD data. From (10), y* is known, the x* is embedded in nonlinear kernel, and thus x* is hard to be solved through this equation. Thus, we need to build the mapping from HD data to LD data and acquire the LD data of new HD data sample fast.We denote HD data samples not containing spatial position information asY~=[y~up1,…,y~upi,…,y~upNup]T∈RNup×D2; then we denote the corresponding LD data as X~=[x~1,…,x~i,…,x~Nup]T∈RNup×q, which can be acquired by GPDM. The mapping can be built through the training model as below:(11)X~T=WDΦ+E,where WD∈Rq×Nk is the weight matrix and Φ∈RNk×Nup is the matrix of Gaussian basis function; thus we have Φk1,k2=ϕ(y~upk2,ck1)=exp⁡(by~upk2-ck12), E∈Rq×Nup is the bias matrix, and Nk≤Nup. b is known; after finishing training of (11), ck1, WD, and E can be confirmed, here, W⌢D and E⌢=B⌢pones(1,Nup) denote the trained values of WD and E, respectively, B⌢p∈Rq×1. Then, we have the mapping of incremental dimension reduction as below:(12)xt*=gyupt*=W⌢DΦyupt*+E⌢,where as to one HD new samples yupt*, let Nup=1, Φ(yupt*)=[ϕ(yupt*,c1),ϕ(yupt*,c2),…,ϕ(yupt*,cNk)]T, and E⌢=B⌢pones(1,1). The mapping as (12) can acquire the corresponding LD data of the new HD data sample; furthermore, it can acquire the multiple corresponding ones of the multiple HD data samples, which only need to add the corresponding column in Φ and E⌢. Equations (5)–(12) constitute IPDRM; however, the learning method of (11) is the key to estimate the human motion in 3D; thus we will discuss it in the next section.In the following, we need to discuss the several advantages of IPDRM. Firstly, the model is used to build the LD space through the known HD data samples, acquire the LD data of the new HD data sample, and provide the human motion global prior information in the LD space. The searching position in the LD space can be confirmed to disambiguate silhouettes through clustering. This is because the mapping from HD space to LD space is one-to-one. Secondly, it can reduce the estimation error and overcome the transient occlusion through linear and nonlinear searching in the LD space via the global prior information; moreover, the HD data samples generated through the large scale nonlinear searching between the two neighbours LD data in the LD space are better and more valid than the large scale nonlinear searching between the two neighbours HD data in the HD space directly, which can be tested as shown in Figure3. Let us see the two experiments (walking and running models); we find that the HD data which are generated through nonlinear searching in the LD space (Figure 3(b)) and converted to human motion models are more like human motion shape than the ones which are generated through the nonlinear searching in the HD space (Figure 3(a)). Finally, searching in the LD space can have lower computation cost than searching in the HD space obviously.Figure 3
Comparison of searching in the HD space and LD space among neighbour data samples.
(a)
Searching in the HD space (some limbs in model are mass and invalid)
(b)
Searching in the LD space (all limbs in model are regular and valid)


## 4.1. Probabilistic Models of GPDM

We give the probabilistic models of the GPDM in [10] as follows:(5)pY∣X,β¯,W=WN2πNDKYDexp⁡-12tr⁡KY-1YW2YT,(6)pX∣α-=px12πN-1qKxqexp⁡-12tr⁡KX-1X2:NX2:NT,In (5), HD data sequence can be denoted by Y=y1,…,yi,…,yNT∈RN×D, yi∈RD, and LD data sequence can be denoted by X=[x1,…,xi,…,xN]T∈RN×q, xi∈Rq. Kernel matrix with parameter β¯=[β1,β2] is KY∈RN×N, KYi,j=kYxi,xj=exp⁡(-β1/2xi-xj2)+β2-1δxixj; W is scale parameter matrix, p(W)=∏m=1D(2/κ2π)exp⁡(-wm2/2κ2), wm>0, κ=10-3; in (6), kernel matrix with parameter α-=[α1,α2,α3,α4] is KX∈RN-1×N-1, (KX)i,j=kX(xi,xj)=α1exp⁡(-α2/2xi-xj2)+α3xiTxj+α4-1δxixj, X1:N-1=[x1,x2,…,xN-1]T, X2:N=[x2,x3,…,xN]T, and x1 conforms Gaussian distribution of q dimensions. α- and β¯ satisfy p(α¯)∝∏iαi-1 and p(β¯)∝∏iβi-1, respectively. Then, we give the conditional distribution as follows:(7)pY*∣X*,Γ=WM2πMDKY*Dexp⁡-12tr⁡KY*-1ZYW2ZYT,(8)pX*∣Γ=px1*2πM-1qKX*qexp⁡-12tr⁡KX*-1ZXZXT,where Γ={X,Y,α-,β-,W}, Y*∈RM×D, X*∈RM×q, ZY=Y*-ATKY-1Y, KY*=B-ATKY-1A, ZX=X2:N*-CTKX-1X2:N, and KX*=D-CTKX-1C. Kernel matrices include A∈RN×M, Ai,j=kY(xi,xj*), B∈RM×M, Bi,j=kY(xi*,xj*), C∈R(N-1)×(M-1), Ci,j=kX(xi,xj*) and D∈R(M-1)×(M-1), and Di,j=kX(xi*,xj*).When training GPDM,Y is known; thus p(Y) is constant, and min⁡(-ln⁡p(X,α¯,β¯,W∣Y))⇔min⁡(-ln⁡p(X,α¯,β¯,W,Y)). Then, the LD data and corresponding parameters can be calculated as below:(9)X,α¯,β¯,W=argmin⁡X,α¯,β¯,W(-ln⁡p(X,α¯,β¯,W,Y))=argmin⁡X,α¯,β¯,W-ln⁡pY∣X,β¯,WpX∣α-pα¯pβ¯pW·pX∣α-pα¯pβ¯pW.Optimization of (9) can use the method of scaled conjugate gradient (SCG) [18, 19]. After the optimization, X can be depicted as shown in Figure 1(e). According to (7), y*~N(μY(x*),σY2(x*)I), the mapping from LD data to HD data can be given by the mean as below:(10)y*=fx*=μY(x*)=YTKY-1kYx1,x*,kYx2,x*,…,kYxN,x*T=YTKY-1kYx*.


## 4.2. The Mapping of Incremental Dimension Reduction

After finishing training, GPDM can not process the new HD data sample, which means it can not acquire the corresponding LD data. From (10), y* is known, the x* is embedded in nonlinear kernel, and thus x* is hard to be solved through this equation. Thus, we need to build the mapping from HD data to LD data and acquire the LD data of new HD data sample fast.We denote HD data samples not containing spatial position information asY~=[y~up1,…,y~upi,…,y~upNup]T∈RNup×D2; then we denote the corresponding LD data as X~=[x~1,…,x~i,…,x~Nup]T∈RNup×q, which can be acquired by GPDM. The mapping can be built through the training model as below:(11)X~T=WDΦ+E,where WD∈Rq×Nk is the weight matrix and Φ∈RNk×Nup is the matrix of Gaussian basis function; thus we have Φk1,k2=ϕ(y~upk2,ck1)=exp⁡(by~upk2-ck12), E∈Rq×Nup is the bias matrix, and Nk≤Nup. b is known; after finishing training of (11), ck1, WD, and E can be confirmed, here, W⌢D and E⌢=B⌢pones(1,Nup) denote the trained values of WD and E, respectively, B⌢p∈Rq×1. Then, we have the mapping of incremental dimension reduction as below:(12)xt*=gyupt*=W⌢DΦyupt*+E⌢,where as to one HD new samples yupt*, let Nup=1, Φ(yupt*)=[ϕ(yupt*,c1),ϕ(yupt*,c2),…,ϕ(yupt*,cNk)]T, and E⌢=B⌢pones(1,1). The mapping as (12) can acquire the corresponding LD data of the new HD data sample; furthermore, it can acquire the multiple corresponding ones of the multiple HD data samples, which only need to add the corresponding column in Φ and E⌢. Equations (5)–(12) constitute IPDRM; however, the learning method of (11) is the key to estimate the human motion in 3D; thus we will discuss it in the next section.In the following, we need to discuss the several advantages of IPDRM. Firstly, the model is used to build the LD space through the known HD data samples, acquire the LD data of the new HD data sample, and provide the human motion global prior information in the LD space. The searching position in the LD space can be confirmed to disambiguate silhouettes through clustering. This is because the mapping from HD space to LD space is one-to-one. Secondly, it can reduce the estimation error and overcome the transient occlusion through linear and nonlinear searching in the LD space via the global prior information; moreover, the HD data samples generated through the large scale nonlinear searching between the two neighbours LD data in the LD space are better and more valid than the large scale nonlinear searching between the two neighbours HD data in the HD space directly, which can be tested as shown in Figure3. Let us see the two experiments (walking and running models); we find that the HD data which are generated through nonlinear searching in the LD space (Figure 3(b)) and converted to human motion models are more like human motion shape than the ones which are generated through the nonlinear searching in the HD space (Figure 3(a)). Finally, searching in the LD space can have lower computation cost than searching in the HD space obviously.Figure 3
Comparison of searching in the HD space and LD space among neighbour data samples.
(a)
Searching in the HD space (some limbs in model are mass and invalid)
(b)
Searching in the LD space (all limbs in model are regular and valid)


## 5. Orthogonal Least Squares Learning of Multiple Outputs

In this section, we propose a method called orthogonal least squares learning of multiple outputs (OLSLMO). In the LD space or HD space, the data are multiple dimensions, for which the feature or model can be described better. Let us denote the model as follows:(13)D=Pθ~+e.Here, matrices of the output and input vectors are, respectively, D=[d1,…,dk,…,dN]T=[d1′,…,dc′,…,dQ1′]∈RN×Q1 and V=[v1,…,vk,…,vN]T∈RN×Q2; P=[p1,…,pk′,…,pM]∈RN×M, Pk,k′=fe(vk,zk′), zk′∈{vk∣k=1,2,…,N}, and θ~=[θ1,…,θk′,…,θM]T=[θ1′,…,θc′,…,θQ1′]∈RM×Q1; error matrix is e=[e1,…,ek,…,eN]T=[e1′,…,ec′,…,eQ1′]∈RN×Q1. Then, the matrix P can be decomposed into this form P=W~A, and W~=[w1,…,wk′,…,wM]∈RN×M is the matrix with orthogonal vector wk′(wk′Twk~′=0, k′≠k~′; wk′Twk~′≠0, k′=k~′; k~′,k′=1,2,…,M), W~TW~∈RM×M is a diagonal matrix; A∈RM×M is the invertible matrix with 1 on the diagonal given by(14)Ai,j=ai,j,i<j,1,i=j,0,i>j;let least squares (LS) estimator θ~=(PTP)-1PTD, g⌢=(W~TW~)-1W~TD, and we can give the derivation as below:(15)Pθ~=W~Aθ~=W~APTP-1PTD=W~AATW~TW~A-1ATW~TD=W~W~TW~-1W~TD=W~g⌢.Equation (13) can be written by(16)D=W~g⌢+e,and let g⌢=[g1,…,gk′,…,gM]T=[g⌢1,…,g⌢c,…,g⌢Q1]∈RM×Q1; according to property of LS estimation, we have(17)min⁡D-Pθ~2N×Q1=min⁡D-W~g⌢2N×Q1=min⁡tr⁡eTeN×Q1=min⁡1Q1∑c=1Q11Nec′Tec′=min⁡⁡1Q1∑c=1Q11Ndc′Tdc′-1N∑k′=1Mwk′Twk′g⌢c2k′,and let Sw={w1,…,wk′,…,wM}, Sw′={w1′,…,wk′′,…,wN′}; minimizing (17) can be achieved by(18)max⁡Sw∑c=1Q1∑k′=1Mwk′Twk′g⌢c2k′dc′Tdc′n~,0<M≤N,Sw∈Sw′,n~=1,2,where gk′=wk′TD/(wk′Twk′); let P′=[p1′,…,pk′′,…,pN′]∈RN×N, Pk,k′′=fe(vk,zk′), and Sw′ is corresponding orthogonal vector set from vector set {p1′,…,pk′′,…,pN′}. The value of M can be confirmed according to the condition D-Pθ~2/(N×Q1)=D-W~g⌢2/(N×Q1)=tr⁡(eTe)/(N×Q1)<ε1, ε1>0. When this condition is satisfied, we let(19)D=Pθ⌢+E⌢′+e′=[P,RI]θ⌢T,BpTT+e′=GF⌢+e′,where RI=ones(N,1)∈RN×1 and error matrix is e′∈RN×Q1. According to the equation and updating method [20] of LS, D⌢ can denote the LS estimator of D, and we let F⌢=(GTG)-1GTD, PG=G(GTG)-1GT, Pθ=P(PTP)-1PT, and P1/θ=I-Pθ and get(20)D⌢=GF⌢=PGD=PθD+P1/θRI(RITP1/θTP1/θRI)-1RITP1/θTD=Pθ~+P1/θRI(RITP1/θTP1/θRI)-1RITP1/θTD,where if condition D-Pθ~2/(N×Q1)=D-W~g⌢2/(N×Q1)=tr⁡(eTe)/(N×Q1)<ε1, ε1>0 is satisfied, and if condition D-GF⌢2/(N×Q1)=D-PGD2/(N×Q1)=tr⁡(e′Te′)/(N×Q1)<ε2<ε1, ε2>0 can also be satisfied; because Pθ is determined, D can get further approximation based on Pθ~ which is playing the leading role during training. The more detailed deviation of the condition can be seen in appendices. Then, we have θ⌢=F⌢1:M,:∈RM×Q1, Bp=F⌢M+1,:∈R1×Q1, and E⌢′=RIBp∈RN×Q1.Then, we summarize the algorithm of OLSLMO as in Algorithm1.Algorithm 1
Input:D, V, n~, b;
Output:zk′,k′=1,2,…,M, θ⌢, Bp;
Fork=1:1:N
Ifk = 1
w
k
s
=
p
s
′, s=1,2,…,N;
Else Ifk ≥ 2
α
i
k
s
=
w
i
T
p
s
′
w
i
T
w
i, s=1,2,…,N,s∉In1:k-1, 1≤i<k;
w
k
s
=
p
s
′
-
∑
i
=
1
k
-
1
α
i
k
s
w
i;
End If
g
k
s
=
(
w
k
s
)
T
D
w
k
s
T
w
k
s;
E
V
s
=
∑
c
=
1
Q
1
w
k
s
T
w
k
s
g
1
s
c
2
d
c
′
T
d
c
′
n
~;
s
*
=
arg
⁡
max
s
⁡
E
V
s; wk=wks*;
I
n
k
=
s
*; zk=vInk;  G=[P′:,In1:k,RI];
F
⌢
=
G
T
G
-
1
G
T
D; θ⌢=F⌢1:k,:;
B
p
=
F
⌢
k
+
1
,
:; E⌢′=RIBp;
IfD-GF⌢2/N×Q1<ε2
Break;
End If
End For
M
=
k;The learning of (11) can use Algorithm 1, we let input parameters D=X~, V=Y~, Q1=q, and Q2=D2 and functions fe(vk,zk′)=ϕ(y~upk2,ck1), vk=y~upk2, k=k2=1,2,…,Nup, and N=Nup, and then get corresponding output parameters ck1=zk′, Nk=M, k1=k′=1,2,…,Nk, W⌢D=θ⌢T, and E⌢=[ones(1,1)Bp]T.


## 6. Low Dimensional Space Incremental Learning

After using low dimensional space searching (LDSS) to estimate some human motion models of corresponding frames in 3D (about one motion cycle), the new HD data samples can be got from the models. The LDSS algorithm is given in the framework of LDSIL shown in Section7.2. Then, continuing to use the mapping of (10) to generate HD data sample will be less accurate, thus incremental learning is needed to update the mapping from LD space to HD space. Equation (10) is so complicated that its parameters can not be updated easily, due to the need of calculation in (9). For this consideration, we still use the model like (11) to build the mapping, which can take advantage of new HD data samples to learn through the method of Algorithm 1. The details of incremental learning will be discussed as follows.The key of incremental learning is how to select the HD data samples to update the mapping from LD space to HD space. The better new samples used to training the mapping can be selected from the estimated HD data with spatial position information removed. The estimated HD data samples contain errors, which can convert into human motion models of corresponding frames. After we get corresponding new HD data samplesY~new*=[y~up1*,…,y~upn**,…,y~upN**]T∈RN*×D2, we can use (12) to obtain the corresponding LD data, (X~new*)T=g(Y~new*T)=[g(y~up1*),…,g(y~upn**),…,g(y~upN**)], and get X~new*=[x~1*,…,x~n**,…,x~N**]T∈RN*×q. We know the LD data is representing the feature of HD data; thus the distance among LD data can describe the similarity of HD data, and we need to select the LD data from X~new*, which are well-distributed in the LD space for getting the corresponding HD data. The selecting rule can be carried out as below. We denote ordinal number as below:(21)fL=i,1≤i≤Nup,bL=i+1,1≤i<Nup,bL=1,i=Nup,and then, we have(22)Dn**=min⁡Dse,Dse∈x~n**-x~fL,x~n**-x~bL,ddddddddddddddddn*=1,2,…,N*,N*>Nup,and the LD data and the corresponding HD data can be selected by(23)Xse*=Xnew*nse,:Yse*=Ynew*nse,:,nse={n*∣Dn**<csxfL-xbL,n*=1,2,…,N*},0<cs≤1,where we let X~se*=[x~s1*,…,x~sis*,…,x~sNs*]T∈RNs×q and Y~se*=[y~s1*,…,y~sis*,…,y~sNs*]T∈RNs×D2 and then get (24)is*=argmin⁡isx~sis*-x~fL-x~sis*-x~bL,is=1,2,…,Ns.Equation (24) means selecting the best well-distributed LD data. Thus, the training data samples of input and output are given by(25)X~neven,:tra=X~bs*=x~1*bs,…,x~ibs*bs,…,x~Nup*bsT∈RNup×q,ix~ibs*bs=x~sis**,Y~neven,:tra=Y~*bs=y~1*bs,…,y~ibs*bs,…,y~Nup*bsT∈RNup×D2,iy~ibs*bs=y~sis**,X~nodd,:tra=X~∈RNup×q,Y~nodd,:tra=Y~∈RNup×D2,ibs=i,nodd=1,3,…,2Nup-1,neven=2,4,…,2Nup,where we let X~tra=[x~1tra,…,x~ntratra,…,x~Nup*tra]T∈RNup*×q and Y~tra=[y~1tra,…,y~ntratra,…,y~Nup*tra]T∈RNup*×D2, X~bs*, Y~bs* are the samples of incremental learning, the X~tra can be seen in the LD space as depicted in Figure 4, and the green ones denote X~bs*; the new mapping can be built through the training model as below: (26)Y~traT=y~1tra,y~2tra,…,y~Nup*tra=y~up1,y~1*bs,y~up2,y~2*bs,…,y~upNup,y~Nup*bs=WD*Φ*+E*=WD*Φ*x~1tra,Φ*x~2tra,…,Φ*x~Nup*tra+E*=WD*x~Nup,Φ*x~Nup*bsΦ*x~1,Φ*x~1*bs,Φ*x~2,Φ*x~2*bs,…,Φ*x~Nup,Φ*x~Nup*bs+E*,where WD*∈RD2×Nk, Φ*∈RNk×Nup*,Φk1,k2*=ϕ*(x~k2tra,ck1*)=exp⁡(bx~k2tra-ck1*2), Φ(xE)=[ϕ(xE,c1*),ϕ(xE,c2*),…,ϕ(xE,cNk*)]T, E*∈RD2×Nup*, Nup*=2Nup, W⌢D* and E⌢*=B⌢p*ones(1,Nup*) denote the trained values of WD* and E*, respectively, and B⌢P*∈RD2×1. After finishing the training of (26), the incremental learning is completed. Then, the mapping is updated, so that the new HD data can be generated through the updated mapping:(27)yupt*=fnewxt*=W⌢D*Φ*xt*+E⌢*,where, similarly, as to the one xt*, let Nup*=1, Φ(xt*)=[ϕ(xt*,c1*),ϕ(xt*,c2*),…,ϕ(xt*,cNk*)]T, and E⌢*=B⌢p*ones(1,1). When we have multiple xt*, we can also get the corresponding multiple yupt*, and it only needs to add the column in Φ* and E⌢*.Figure 4
The low dimensional data through increment learning.Training equation (26) can use Algorithm 1; we let input parameters D=Y~tra, V=X~tra, Q1=D2, and Q2=q and functions fe(vk,zk′)=ϕ*(x~k2tra,ck1*), vk=x~upk2, k=k2=1,2,…,Nup*, and N=Nup* and then get corresponding output parameters ck1*=zk′, Nk=M, k1=k′=1,2,…,Nk, W⌢D*=θ⌢T, and E⌢*=[ones(1,1)Bp]T. Let us summarize the incremental learning algorithm  as in Algorithm 2.Algorithm 2
Input:X~new*, Y~new*, b, cs;
Output:ck1*,k1=1,2,…,Nk, W⌢D*, E⌢*, X~tra;
Fori=1:1:Nup
Ifi = Nup
f
L
=
i; bL=1;
Else
f
L
=
i; bL=i+1;
End If
k
D
=
1;
Forn*=1:1:N*
D
n
*
*
=
min
⁡
D
se
,
D
se
∈
x
~
n
*
*
-
x
~
f
L
,
x
~
n
*
*
-
x
~
b
L;
IfDn**<csxfL-xbL
I
N
k
D
=
n
*;
End If
k
D
=
k
D
+
1;
End For
X
se
*
=
X
new
*
I
N
1
:
k
D
,
:
;
Y
se
*
=
Y
new
*
I
N
1
:
k
D
,
:;
Foris=1:1:Ns
D
gap
⁡
i
s
=


s
x
~
i
s
*
-
x
~
f
L
-


s
x
~
i
s
*
-
x
~
b
L;
End For
i
s
*
=
arg
⁡
min
i
s
⁡
D
gap
⁡
i
s
;
i
bs
=
i; bsx~ibs*=sx~is**;bsy~ibs*=sy~is**;
End For
X
~
n
even
,
:
tra
= bsX~; Y~neven,:tra=bsY~*; X~nodd,:tra=X~; Y~nodd,:tra=Y~;
D
=
Y
~
tra; V=X~tra; fevk,zk′=ϕ*x~k2tra,ck1*;
Q
1
=
D
2; Q2=q; vk=x~k2tra; zk′=ck1*; k=k2; k′=k1;
Use Algorithm1 to get zk′,k′=1,2,…,M,θ⌢ and Bp;
c
k
1
*
=
z
k
′, k1=k′=1,2,…,Nk, Nk=M; W⌢D*=θ⌢T;
E
⌢
*
=
ones
1,1
B
p
T;The mapping of (27) is incremental learn mapping, which replaces (10) of IPDRM; it can promote the second advantage of IPDRM mentioned above, due to the mapping from LD data to HD data being updated to get more accurate. Worthy of attention is that at least one x~n** must distribute between x~i and x~i+1, so that effectiveness of incremental learning can be better, thus using LDSS to generate HD data samples should follow this.


## 7. Human Motion Estimation via Searching in Low Dimensional Space

### 7.1. The Method of Searching

We estimate the human motion model via searching in low dimensional space. The method of searching needs to combine linear searching and nonlinear searching in order to generate the better corresponding HD data sample, which is depicted as Figure5. In Figure 5, the hollow dots denote the new LD data acquired by linear searching and nonlinear searching in the LD space; the solid dots denote the known LD data acquired by dimension reduction of known HD data samples. Let us discuss more details as follow.Figure 5
The method of searching.We takeX~′=[x~1′,…,x~i′,…,x~Nup′′]T∈RNup′×q (X~′=X~ or X~tra), Y~′=[y~1′,…,y~i′,…,y~Nup′′]T∈RNup′×D2 (Y~′=Y~ or Y~tra), and the mapping fs (fs=f or fnew), for example, to discuss the method. As to the new HD data sample yupt*, we use mapping g to obtain the xt* (xt*=g(yupt*)). Then, the clustering of the xt* can be executed through the comparison of distance to find the nearest known LD data:(28)i*=argmin⁡ixt*-x~i′,i=1,2,…,Nup′.After finding the x~i*′, Nls denotes the amount of new searched LD data; we need linear searching to get the new LD data between the two known neighbours LD data near x~i*′, like this(29)sl1=x~B2′+l(x~B1′-x~B2′)(Nls+1),sl2=x~B1′+l(x~i*′-x~B1′)Nls+1,sl3=x~i*′+l(x~F1′-x~i*′)(Nls+1),sl4=x~F1′+l(x~F2′-x~F1′)(Nls+1),dddddddl=1,2,…,Nls;here, Fns=i*+ns, Bns=i*-ns, and ns∈N+; let ns=1,2; then we have Fns=Fns-Nup′, Fns>Nup′, and Bns=Bns+Nup′, Bns<1. Besides, we also need nonlinear searching to get them similarly, and Nlns denotes the number of new ones(30)s⌢l′1=x~B1′+u~l′1,s⌢l′2=x~i*′+u~l′2,s⌢l′3=x~F1′+u~l′3,l′=1,2,…,Nlns,where u~l′sv=η~l′L~(u~l′sv/u~l′sv), sv=1,2,3,…, L~=argmaxD~, D~∈{x~i*′-x~B1′,x~B1′-x~B2′,x~i*′-x~F1′,x~F1′-x~F2′},  η~l′=rand(1,1), and u~l′sv=rand(q,1). From (29) and (30), the set containing searched LD data can be(31)SL=x~B2′,s11,…,sNls1,x~B1′,s12,…,sNls2,x~i*′,s13,…,sNls3,s⌢Nlns3x~F1′,s14,…,sNls4,x~F2′,s⌢11,…,s⌢Nlns1,s⌢12,…,s⌢Nlns2,s⌢13,…,s⌢Nlns3;use the mapping fs to get another set as below:(32)SLH⁡=yupt*,ya′=fsxa∣xa∈SL;according to (2), we have (33)y⌢a*=argmax⁡y⌢apmt∣y⌢p*,y⌢a,y⌢a∈SLH⁡,where the optimal prior information of the pose y⌢a*∈RD2 can be obtained, and then the optimal prior HD data can be ypri=[y⌢p*T,y⌢a*T]T∈RD1, and y⌢p* is the optimal spatial position information which can be obtained by SEMAS.On the basis of work above, we can estimate the human motion model from multiview silhouettes through Bayesian theory. Whenypri is known, p(ypti∣mt) is unknown and Ne′ is large enough; we let p(ypti)=q(ypti∣y-=ypri), and the estimation can be achieved by(34)wbi=p(ypti,mt)qypti∣y-=ypri=pmt∣yptipyptiqypti∣y-=ypri=pmt∣ypti=wti,(35)ypt*=Eypti∣mt=∫yptipypti∣mtdypti≈∑i=1Ne′yptipypti∣mt=∑i=1Ne′yptipypti,mtp(mt)=∑i=1Ne′yptipypti,mt∑i=1Ne′pypti,mt=∑i=1Ne′yptiwtiqypti∣y-=ypri∑i=1Ne′wtiqypti∣y-=ypri≈1/Ne′∑i=1Ne′yptiwti1/Ne′∑i=1Ne′wti=∑i=1Ne′w~tiypti,and, after getting ypt*, the human motion model can be drawn in the space. Worthy of attention is that q(ypti∣ypri) is multivariate normal distribution density function, which takes y-=ypri as the mean, Ne′ is the number of generating ypti, ypri is contained as one of Ne′ samples, which can let yptNe′=ypri, and w~ti=wti/∑i=1Ne′wti is normalized weight. The method of searching before or after incremental learning is also like this; it only needs to replace some variables to be carried out.The poses of human motion are continuous, which are shown by the human motion models. The pose estimated last time is close to the one of this time, thus the new HD sampleyupt* can be seen as the pose estimated last time, and the corresponding LD data xt* can be acquired by IPDRM. Then, we find the other LD data close to xt* through clustering to confirm the searching position in the LD space, so that the HD data samples closing to the true data this time can be generated efficiently through searching in the nearby area of this position. For these advantages, the proposal algorithm (LDSIL) can disambiguate silhouettes, overcome the transient occlusion, and reduce estimation error. Furthermore, the performance can be promoted through its incremental learning which is used to update the mapping from the LD space to the HD space.As to the estimation of Bayesian theory above, we use SEMAS and the proposed method of searching to find the optimal prior HD dataypri, which can make p(ypti,mt)/q(ypti∣y-)=p(mt|ypti) large enough. When ypti=ypri, according to (34), somehow, p(ypti,mt) also become large enough. Then, it means that ypri can be seen as the sample generating from the unimodal distribution density p(ypti∣mt), due to p(ypti,mt)∝p(ypti∣mt)=p(ypti,mt)/p(mt), and p(mt) is known. Moreover, the true HD data can be seen as the mean relative to the distribution density p(ypti∣mt), and thus its value of distribution density will be larger than others. On the whole, it also means that ypri is the sample which is close to the true HD data; thus taking ypri as the mean (y-=ypri) with the small scale variance to search the HD data in the HD space for only once can generate the HD data samples which will also have the large value of distribution density p(ypti∣mt); then the mean of these samples can be close to mean relative to distribution density p(ypti∣mt); in other words, it will be close to the true HD data. The work above can avoid the small scale searching (sampling) the HD data for many times or the large scale searching (sampling) the HD data directly and generating the invalid HD data to estimate the human motion in 3D.Then, we will discuss more details about the whole procedure of human motion estimation in the following section.


### 7.2. The Procedure of Algorithm for Human Motion Estimation

The work from Sections2–7 can be summarized as the complete algorithm procedure of LDSIL to estimate human motion in 3D. The framework of LDSIL can be seen as shown in Figure 6. Let us give more detailed description for the framework.Figure 6
Framework of LDSIL.From Figure6, we can see the framework of LDSIL. Firstly, LDSIL is using the LDSS to estimate the human motion in 3D based on the small amount of HD data samples and get the estimated HD data covering about one motion cycle. The framework of LDSS consists of SEMAS, linear and nonlinear searching in LD space, IPDRM, and Bayesian theory. The estimated HD data with the spatial position information removed can be selected as the new HD data samples and used for updating the mapping from LD space to HD space through IPDRM. Secondly, the IPDRM is used to obtain the corresponding LD data of the new HD data samples through its mapping g of incremental dimension reduction. Selecting the new HD data samples can be achieved through comparing the distance among the corresponding LD data. After the comparison is finished, the LD data and the corresponding HD data can be obtained and used to update the mapping from LD space to HD space for the incremental learning of LD space. Thirdly, the new mapping fnew updated by the selected HD data and the corresponding LD data can generate the HD data more accurately; thus the estimation can achieve a better performance with the help of new mapping fnew. Finally, the estimation of human motion in 3D from the multiview silhouettes is carried out through the framework including SEMAS, linear and nonlinear searching in LD space, IPDRM with new mapping fnew, and Bayesian theory. However, when LDSIL is used to estimate the human motion, these data are known as follows.We have the sequence of the small amount of HD data samplesY~ (no spatial position information) covering one motion cycle and the HD data of initial frame (1st frame) containing the spatial position information (6 dimensions in total) yp0*=yp1. Ft can denote the number of estimatedframes without incremental learning, which are used for the sample selection of incremental learning, and Ft* can denote the number of estimated frames with incremental learning. Then, we begin the procedure as follows.LDSIL for Human Motion Estimation
(1)
Use (9) to reduce the dimension of Y~, obtain X~, build the mapping f, and let yp0*=yp1, t=1.
(2)
Use Algorithm1 to training equation (11), build the mapping g, get IPDRM, and let mapping fs=f.
(3)
When estimating thetth frame, use SEMAS and set relevant parameters to find the spatial position information. Let y0pk⌢1=ypt-1*(1:2), select u~ti=ηiL(uti/uti), and get the optimal planimetric position information y1p*=yp*∈R2. Then, let y0pk⌢1=ypt-1*(3:6), select u~ti=Luti, and get the optimal height and rotation angle information y2p*=yp*∈R4. At last, obtain y⌢p*=[y1p*T,y2p*T]T.
(4)
Letyupt*=ypt-1*(7:D1), obtain the LD data of the tth frame xt*=g(yupt*), and calculate i*=argmin⁡i⁡xt*-x~i′,  i=1,2,…,Nup. Search the data in the LD space, get the set SLH⁡={yupt*,ya′=fsxa∣xa∈SL}, and calculate y⌢a*=argmax⁡y⌢a⁡p(mt∣y⌢p*,y⌢a), y⌢a∈SLH⁡, and then get ypri=[y⌢p*T,y⌢a*T]T.
(5)
Generatepyti, i=1,2,…,Ne′, through q(ypti∣y-=ypri), let yptNe′=ypri, calculate w~ti=wti/∑i=1Ne′wti, i=1,2,…,Ne′, and get ypt*=∑i=1Ne′w~tiypti. Then, human motion model can be outputted through ypt*, and let y~upn**=ypt*(7:D1), n*=t. If t≤Ft, t=t+1, return to step (3); otherwise, get the estimated HD data Y~new* (new HD data samples) and go to next step.
(6)
If needing incremental learning, go to next step; otherwise end this algorithm.
(7)
After gettingY~new*, use mapping g to obtain (X~new*)T=g(Y~new*T)=[g(y~up1*),…,g(y~upn**),…,g(y~upN**)]. Then, get the training samples Y~tra and X~tra selected from Y~new* and X~new* by the method proposed above, and use Algorithm 2 to carry out incremental learning to train (26). After finishing the training, get updated mapping fnew. Let X~′=X~tra, Ft=Ft*, yp0*=yp1, and mapping fs=fnew, and return to step (3).


## 7.1. The Method of Searching

We estimate the human motion model via searching in low dimensional space. The method of searching needs to combine linear searching and nonlinear searching in order to generate the better corresponding HD data sample, which is depicted as Figure5. In Figure 5, the hollow dots denote the new LD data acquired by linear searching and nonlinear searching in the LD space; the solid dots denote the known LD data acquired by dimension reduction of known HD data samples. Let us discuss more details as follow.Figure 5
The method of searching.We takeX~′=[x~1′,…,x~i′,…,x~Nup′′]T∈RNup′×q (X~′=X~ or X~tra), Y~′=[y~1′,…,y~i′,…,y~Nup′′]T∈RNup′×D2 (Y~′=Y~ or Y~tra), and the mapping fs (fs=f or fnew), for example, to discuss the method. As to the new HD data sample yupt*, we use mapping g to obtain the xt* (xt*=g(yupt*)). Then, the clustering of the xt* can be executed through the comparison of distance to find the nearest known LD data:(28)i*=argmin⁡ixt*-x~i′,i=1,2,…,Nup′.After finding the x~i*′, Nls denotes the amount of new searched LD data; we need linear searching to get the new LD data between the two known neighbours LD data near x~i*′, like this(29)sl1=x~B2′+l(x~B1′-x~B2′)(Nls+1),sl2=x~B1′+l(x~i*′-x~B1′)Nls+1,sl3=x~i*′+l(x~F1′-x~i*′)(Nls+1),sl4=x~F1′+l(x~F2′-x~F1′)(Nls+1),dddddddl=1,2,…,Nls;here, Fns=i*+ns, Bns=i*-ns, and ns∈N+; let ns=1,2; then we have Fns=Fns-Nup′, Fns>Nup′, and Bns=Bns+Nup′, Bns<1. Besides, we also need nonlinear searching to get them similarly, and Nlns denotes the number of new ones(30)s⌢l′1=x~B1′+u~l′1,s⌢l′2=x~i*′+u~l′2,s⌢l′3=x~F1′+u~l′3,l′=1,2,…,Nlns,where u~l′sv=η~l′L~(u~l′sv/u~l′sv), sv=1,2,3,…, L~=argmaxD~, D~∈{x~i*′-x~B1′,x~B1′-x~B2′,x~i*′-x~F1′,x~F1′-x~F2′},  η~l′=rand(1,1), and u~l′sv=rand(q,1). From (29) and (30), the set containing searched LD data can be(31)SL=x~B2′,s11,…,sNls1,x~B1′,s12,…,sNls2,x~i*′,s13,…,sNls3,s⌢Nlns3x~F1′,s14,…,sNls4,x~F2′,s⌢11,…,s⌢Nlns1,s⌢12,…,s⌢Nlns2,s⌢13,…,s⌢Nlns3;use the mapping fs to get another set as below:(32)SLH⁡=yupt*,ya′=fsxa∣xa∈SL;according to (2), we have (33)y⌢a*=argmax⁡y⌢apmt∣y⌢p*,y⌢a,y⌢a∈SLH⁡,where the optimal prior information of the pose y⌢a*∈RD2 can be obtained, and then the optimal prior HD data can be ypri=[y⌢p*T,y⌢a*T]T∈RD1, and y⌢p* is the optimal spatial position information which can be obtained by SEMAS.On the basis of work above, we can estimate the human motion model from multiview silhouettes through Bayesian theory. Whenypri is known, p(ypti∣mt) is unknown and Ne′ is large enough; we let p(ypti)=q(ypti∣y-=ypri), and the estimation can be achieved by(34)wbi=p(ypti,mt)qypti∣y-=ypri=pmt∣yptipyptiqypti∣y-=ypri=pmt∣ypti=wti,(35)ypt*=Eypti∣mt=∫yptipypti∣mtdypti≈∑i=1Ne′yptipypti∣mt=∑i=1Ne′yptipypti,mtp(mt)=∑i=1Ne′yptipypti,mt∑i=1Ne′pypti,mt=∑i=1Ne′yptiwtiqypti∣y-=ypri∑i=1Ne′wtiqypti∣y-=ypri≈1/Ne′∑i=1Ne′yptiwti1/Ne′∑i=1Ne′wti=∑i=1Ne′w~tiypti,and, after getting ypt*, the human motion model can be drawn in the space. Worthy of attention is that q(ypti∣ypri) is multivariate normal distribution density function, which takes y-=ypri as the mean, Ne′ is the number of generating ypti, ypri is contained as one of Ne′ samples, which can let yptNe′=ypri, and w~ti=wti/∑i=1Ne′wti is normalized weight. The method of searching before or after incremental learning is also like this; it only needs to replace some variables to be carried out.The poses of human motion are continuous, which are shown by the human motion models. The pose estimated last time is close to the one of this time, thus the new HD sampleyupt* can be seen as the pose estimated last time, and the corresponding LD data xt* can be acquired by IPDRM. Then, we find the other LD data close to xt* through clustering to confirm the searching position in the LD space, so that the HD data samples closing to the true data this time can be generated efficiently through searching in the nearby area of this position. For these advantages, the proposal algorithm (LDSIL) can disambiguate silhouettes, overcome the transient occlusion, and reduce estimation error. Furthermore, the performance can be promoted through its incremental learning which is used to update the mapping from the LD space to the HD space.As to the estimation of Bayesian theory above, we use SEMAS and the proposed method of searching to find the optimal prior HD dataypri, which can make p(ypti,mt)/q(ypti∣y-)=p(mt|ypti) large enough. When ypti=ypri, according to (34), somehow, p(ypti,mt) also become large enough. Then, it means that ypri can be seen as the sample generating from the unimodal distribution density p(ypti∣mt), due to p(ypti,mt)∝p(ypti∣mt)=p(ypti,mt)/p(mt), and p(mt) is known. Moreover, the true HD data can be seen as the mean relative to the distribution density p(ypti∣mt), and thus its value of distribution density will be larger than others. On the whole, it also means that ypri is the sample which is close to the true HD data; thus taking ypri as the mean (y-=ypri) with the small scale variance to search the HD data in the HD space for only once can generate the HD data samples which will also have the large value of distribution density p(ypti∣mt); then the mean of these samples can be close to mean relative to distribution density p(ypti∣mt); in other words, it will be close to the true HD data. The work above can avoid the small scale searching (sampling) the HD data for many times or the large scale searching (sampling) the HD data directly and generating the invalid HD data to estimate the human motion in 3D.Then, we will discuss more details about the whole procedure of human motion estimation in the following section.


## 7.2. The Procedure of Algorithm for Human Motion Estimation

The work from Sections2–7 can be summarized as the complete algorithm procedure of LDSIL to estimate human motion in 3D. The framework of LDSIL can be seen as shown in Figure 6. Let us give more detailed description for the framework.Figure 6
Framework of LDSIL.From Figure6, we can see the framework of LDSIL. Firstly, LDSIL is using the LDSS to estimate the human motion in 3D based on the small amount of HD data samples and get the estimated HD data covering about one motion cycle. The framework of LDSS consists of SEMAS, linear and nonlinear searching in LD space, IPDRM, and Bayesian theory. The estimated HD data with the spatial position information removed can be selected as the new HD data samples and used for updating the mapping from LD space to HD space through IPDRM. Secondly, the IPDRM is used to obtain the corresponding LD data of the new HD data samples through its mapping g of incremental dimension reduction. Selecting the new HD data samples can be achieved through comparing the distance among the corresponding LD data. After the comparison is finished, the LD data and the corresponding HD data can be obtained and used to update the mapping from LD space to HD space for the incremental learning of LD space. Thirdly, the new mapping fnew updated by the selected HD data and the corresponding LD data can generate the HD data more accurately; thus the estimation can achieve a better performance with the help of new mapping fnew. Finally, the estimation of human motion in 3D from the multiview silhouettes is carried out through the framework including SEMAS, linear and nonlinear searching in LD space, IPDRM with new mapping fnew, and Bayesian theory. However, when LDSIL is used to estimate the human motion, these data are known as follows.We have the sequence of the small amount of HD data samplesY~ (no spatial position information) covering one motion cycle and the HD data of initial frame (1st frame) containing the spatial position information (6 dimensions in total) yp0*=yp1. Ft can denote the number of estimatedframes without incremental learning, which are used for the sample selection of incremental learning, and Ft* can denote the number of estimated frames with incremental learning. Then, we begin the procedure as follows.LDSIL for Human Motion Estimation
(1)
Use (9) to reduce the dimension of Y~, obtain X~, build the mapping f, and let yp0*=yp1, t=1.
(2)
Use Algorithm1 to training equation (11), build the mapping g, get IPDRM, and let mapping fs=f.
(3)
When estimating thetth frame, use SEMAS and set relevant parameters to find the spatial position information. Let y0pk⌢1=ypt-1*(1:2), select u~ti=ηiL(uti/uti), and get the optimal planimetric position information y1p*=yp*∈R2. Then, let y0pk⌢1=ypt-1*(3:6), select u~ti=Luti, and get the optimal height and rotation angle information y2p*=yp*∈R4. At last, obtain y⌢p*=[y1p*T,y2p*T]T.
(4)
Letyupt*=ypt-1*(7:D1), obtain the LD data of the tth frame xt*=g(yupt*), and calculate i*=argmin⁡i⁡xt*-x~i′,  i=1,2,…,Nup. Search the data in the LD space, get the set SLH⁡={yupt*,ya′=fsxa∣xa∈SL}, and calculate y⌢a*=argmax⁡y⌢a⁡p(mt∣y⌢p*,y⌢a), y⌢a∈SLH⁡, and then get ypri=[y⌢p*T,y⌢a*T]T.
(5)
Generatepyti, i=1,2,…,Ne′, through q(ypti∣y-=ypri), let yptNe′=ypri, calculate w~ti=wti/∑i=1Ne′wti, i=1,2,…,Ne′, and get ypt*=∑i=1Ne′w~tiypti. Then, human motion model can be outputted through ypt*, and let y~upn**=ypt*(7:D1), n*=t. If t≤Ft, t=t+1, return to step (3); otherwise, get the estimated HD data Y~new* (new HD data samples) and go to next step.
(6)
If needing incremental learning, go to next step; otherwise end this algorithm.
(7)
After gettingY~new*, use mapping g to obtain (X~new*)T=g(Y~new*T)=[g(y~up1*),…,g(y~upn**),…,g(y~upN**)]. Then, get the training samples Y~tra and X~tra selected from Y~new* and X~new* by the method proposed above, and use Algorithm 2 to carry out incremental learning to train (26). After finishing the training, get updated mapping fnew. Let X~′=X~tra, Ft=Ft*, yp0*=yp1, and mapping fs=fnew, and return to step (3).


## 8. Experiments and Evaluations

We tested the performance of our proposed algorithm (LDSIL) in three views (Nv=3) by comparison with other three traditional algorithms, which included annealed particle filter [4] (APF), Gaussian particle filter [21] (GPF), and particle filter [22] (PF). The tested performance included disambiguating silhouettes, overcoming the transient occlusion, and estimation error. Besides, we compared LDSIL with LDSS to test the validity of incremental learning independently, and the LDSS was similar to LDSIL somehow, so it had the performance of LDSIL, but it was without the incremental learning, and the mapping from LD space to HD space was not updated, which still used (10) to generate HD data. The performance of LDSIL and LDSS could be seen by the comparison of the mean error and maximum error.Firstly, we were testing disambiguating silhouettes and taking walking and running motion as the test cases. Then, we saw in the subfigures of Figure7, after estimating several frames, APF, GPF, and PF could not disambiguate silhouettes to estimate the human motion in 3D; on the contrary, LDSIL could disambiguate silhouettes to estimate the one that was close to the true data. The ambiguity meant a set of silhouettes such that the pictures of Figure 1(c) could not distinguish positions of the limbs from the human body; for example, the silhouettes could not depict which one of the legs was in front or in back. Thus, we could see in the subfigures of Figures 7(a) and 7(b) that the estimated human motion projections in all views from APF, GPF, and PF showed that the positions of left leg (white line) and right leg (gray line) were, respectively, opposite by the comparison of true data. The corresponding human motion models were reflecting the same phenomenon, and the left leg (gray one) and right leg (black one) of the model were also opposite.Figure 7
Several algorithms were disambiguating silhouettes to estimate the human motion in 3D.
(a)
Walking
(b)
RunningSecondly, we tested overcoming transient occlusion. The transient occlusion was equal to estimate the adjacent and discontinuous frame, due to transient occlusion of cameras from all the views, which would lead to the result that the silhouettes from all the views could not be obtained. The two frames (initial data and true data) between which the interval of 10 frames existed were used to test all the algorithms. We could find in the subfigures of Figures8(a) and 8(b) that APF, GPF, and PF could not estimate the human motion models accurately, which showed that the limbs from the human body mismatched the images of all the views. However, LDSIL could estimate the models more accurate than APF, GPF, and PF, which were depicted in the corresponding subfigures, and the limbs from the human body matched the images of all the views better like the true data.Figure 8
Several algorithms were overcoming transient occlusion to estimate the human motion in 3D.
(a)
Walking
(b)
RunningThirdly, we used 50 continuous frames to test the estimation error of each algorithm. We used walking and running motion sequences (walking 1–3, running 1-2) to test, whose spatial positions had obvious change. The subfigures of Figure9 showed the experiment results, the mean error and standard deviation from LDSIL were the smallest among these algorithms, and the maximum error shown in the identifications of the corresponding subfigures was also the smallest; besides, the errors of most frames from LDSIL were smallest. The error of each frame could be computed as below [23]: (36)Er=1Mar∑ma=1MarBxma-Bxma′,Bxma,Bxma′∈R3,where Bxma and Bxma′ were, respectively, true position and estimated position of the joint marker in the model and Mar was the number of the markers. From Figures 9(b)–9(f), we found that the errors of all the tested algorithms were close in the estimation of the first 20 frames, but the error from LDSIL could be kept smaller than APF, GPF, and PF after the 20th frame. However, it was reasonable that the errors from APF, GPF, and PF might be close to LDSIL in some frame estimation, because initial data was close to the true data and the probability was stochastic, which were still larger than the LDSIL on the whole.Figure 9
The estimation errors of several algorithms.
(a)
Mean error and standard deviation
(b)
The error and maximum error in walking 1
(c)
The error and maximum error in running 1
(d)
The error and maximum error in walking 2
(e)
The error and maximum error in running 2
(f)
The error and maximum error in walking 3Finally, we tested the validity of incremental learning of LDSIL. Compared with LDSIL, LDSS had no incremental learning. We selected 3 motion sequences to test the estimation error for 6 times arbitrarily, which included known and unknown spatial position information. We were setting the same parameters in the experiments, and the performance of LDSIL was better than LDSS as shown in the Table1, according to the comparison of the mean errors and maximum error. The experiment results indicated that the incremental learning played the key role, because LDSIL had the smaller mean error and maximum error in the tests.Table 1
The test of incremental learning.
Algorithm
Sequence (total: 30 frames)
Mean error (mm)
Maximum error (mm)
Spatial position information
LDSS
1
54.87
77.24
Unknown
2
60.89
109.90
Unknown
3
34.23
44.02
Unknown
1
54.01
146.93
Known
2
28.08
39.63
Known
3
28.23
37.79
Known
LDSIL
1
49.11
73.47
Unknown
2
43.03
74.06
Unknown
3
31.52
39.60
Unknown
1
25.54
43.42
Known
2
24.36
35.74
Known
3
21.24
29.84
Known


## 9. Conclusions

We could see the experiments above, which indicated that the proposed algorithm (LDSIL) could contribute a lot to estimate the human motion containing the spatial position information in 3D from the multiview silhouettes. The results showed that LDSIL had the best performance including disambiguating silhouettes, overcoming the transient occlusion, and reducing estimation error by the comparison of the other three algorithms (APF, GPF, and PF). Meanwhile, the feasibility and performance of incremental learning in the LDSIL were also validated by the experiments. In addition, the segmentation of multiview images must achieve the high quality; otherwise the results in the experiments would be affected. Our work in this paper had some limitation as follows. Firstly, the initialized frame of the human motion model (the 1st frame) was known, which contained the spatial position information; secondly, the spatial position information and pose of the estimated human motion in 3D were changing regularly, and the human motion model had markers. In the future, these limitations would be improved. Our work would solve the problem that initialized frame needed to be known and focus on complex human motion estimation. Moreover, the human motion model would have no marker [24, 25].



---

*Source: 671419-2015-03-16.xml*