AAM and Non-rigid Registration in Augmented Reality

1. Introduction

Augmented Reality copes with the problem of dynamically augmenting or enhancing the real world with computer generated virtual objects [Azuma, 1997; Azuma, 2001]. Registration is one of the most pivotal problems in augmented reality applications. Typical augmented reality applications track 2D patterns on rigid planar objects in order to acquire the pose of the camera in the scene. Although the problem of rigid registration has been widely studied [Yuan et al., 2005; Yuan et al., 2006; Guan et al., 2008a; Guan et al., 2008b; Li et al., 2008], non-rigid registration is recently receiving more and more attention. There are many non-rigid objects existing in the real world such as animated faces, deformable cloth, hand and so forth. How to overlaid virtual objects on the non-rigid objects is particular challenging.

Recently, many related non-rigid registration approaches have been reported. In many cases (e.g. human faces), only a few feature points can be reliably tracked. In [Bartoli et al., 2004], a non-rigid registration method using point and curve correspondences was proposed to solve this problem. They introduced curves into the non-rigid factorization framework because there are several curves (e.g. the hairline, the eyebows) can be used to determine the mapping. The mapping function is computed from both point and curve correspondences. This method can successfully augment the non-rigid object with a virtual object. In [Pilet et al., 2005; Pilet et al., 2007], they presented a real-time method for detecting deformable surfaces with no need whatsoever for a prior pose knowledge. The deformable 2D meshes are introduced. With the use of fast wide baseline matching algorithm, they can superimpose an appropriately deformed logo on the T-shirt. These methods are robust to large deformations, lighting changes, motion blur and occlusions. To align the virtual objects generated by computers with the real world seamlessly, the accurate registration data should be provided. In general, registration can be achieved by solving a point matching problem. The problem is to find the correspondence between two sets of tracked feature points. Therefore, the detection of feature points and the points tracking are the two main problems. Rigid object detection and tracking have been extensively studied and effective, robust, and real-time solutions proposed [Lowe, 2004; Lepetit & Fua, 2005; Lepetit et al., 2005; Rosten & Drummond, 2005]. Non-rigid object detection and tracking is far more complex because the object is deformable and not only the registration data but also a large number of deformation parameters must be estimated.

Active Appearance Models (AAM), introduced a few years ago [Cootes et al., 1998; Cootes et al., 2001], are commonly used to track non-rigid objects, such as faces and hands. There are many methods have been proposed to track non-rigid objects using AAM. A working system for finding and tracking a human face and its features using active appearance models was presented in [Ahlberg, 2001]. A wireframe model is adapted to the face in each image. Then the model is matched to the face in the input image using active appearance algorithm. In [Sung & Kim, 2004], the previous 2D AAM is extended to 3D shape model and modified model fitting algorithm was proposed. In [Markin & Prakash, 2006], occluded images are included into AAM training data to solve the occlusion and self-occlusion problem. This approach can improve the fitting quality of the algorithm.

With known coordinates of 3D points in the world coordinates and the corresponding 2D image points, the camera pose can be estimated. For non-rigid objects, AAM algorithm is a robust method to acquire the 2D image points. It has been proven to be a useful method for matching any of the statistical models to a new image rapidly. The 3D points of the non-rigid objects can be represented by a linear combination of a set of 3D basis shapes. By varying the configuration weights and camera pose parameters, the error between the estimated 2D points (projected by the estimated 3D shapes using estimated camera pose parameters) and 2D tracking points can be minimized.

Many methods have been proposed to recover 3D shape basis from 2D image sequences. In [Tomasi & Kanade, 1992], the factorization method is used to recover shape and motion from a sequence of images under orthographic projection. The image sequence is represented as a measurement matrix. It is proved that under orthography, the measurement matrix is of rank 3 and can be factored into 3D pose and 3D shape matrix. Unfortunately this technique can not be applied to non-rigid deforming objects, since they are based on the rigidity assumption. The technique based on a non-rigid model is proposed to recover 3D non-rigid shape models under scaled orthographic projection [Bregler et al., 2000]. The 3D shape in each frame can be expressed by a linear combination of a set of K basis shapes. Under this model, the 2D tracking matrix is of rank 3K and can be factored into 3D pose, object configuration and 3D basis shapes with the use of SVD.

In this chapter, a novel non-rigid registration method for augmented reality applications with the use of AAM and factorization method is introduced. We focus on AAM algorithm and factorization method which can obtain the 3D shape basis, object configuration and 3D pose simultaneously. The following demonstrations are mainly based on the researches presented in [Tomasi & Kanade, 1992; Cootes et al., 1998; Bregler et al., 2000; Cootes et al., 2001; Xiao et al., 2004; Zhu et al., 2006; Tian et al., 2008]. In section 2, we will have a detailed review of active appearance models and focus on the way of how to track non-rigid objects. In section 3, we will illustrate how to compute the 3D shape basis, the camera rotation matrix and configuration weights of each training frame simultaneously from the 2D tracking data using factorization algorithm. In section 4, we will introduce how to compute the precise configuration weights and camera rotation matrix by optimization method and the experimental results are also given.

2. Tracking non-rigid objects using AAM

Tracking non-rigid objects using AAM includes five main steps. The first step is to obtain landmarks in training image set. Then establish the shape model and texture model separately. These two models are unified into one appearance model in the next step. Finally the difference between the closest synthesized AAM model and input image is minimized to get the accurate tracking result. The flowchart is shown in Figure 1.

Figure 1.

The flowchart of tracking non-rigid objects using AAM.

2.1. Obtain landmarks

Figure 2.

Examples of traning images manually labeled with consistent landmarks.

Before establishing the models, we should place hundreds of points in each 2D training image. The landmarks in each image should be consistent. We usually choose the points of high curve or junctions as landmark points which will control the shape of the target object strictly. The acquisition process is cumbersome especially when it is done manually. Researchers have looked for different ways to reduce the burden. Ideally, one only need to place points in one image and the corresponding points in the remain traning images can be found automatically. Obviously, this is impossible. However many semi-automatic methods have been proposed and can successfully find the correspondences accross an image set. Here we focus on the tracking procedure, more details about semi-autometic placement of landmarks can be found in [Sclaroff & Pentland, 1995; Duta et al., 1999; Walker et al., 2000; Andresen & Nielsen, 2001]. The examples of training images manually labeled with consistent landmarks are shown in Figure 2. We use 7 training images taken from different viewpoints. Each image is labeled with 36 landmarks shown as the “ $$\times$$ ” in the figures.

Given a training image set with key landmark points are marked on each example object, we can establish the shape and texture variations. These approaches will be detailedly illustrated in the following sections.

2.2. Establish shape model

Suppose L denotes the number of shapes in a training image set and m is the number of key landmark points on each image. The vector representation for each shape would then be:

$$X_i={\left[x_1,x_2,\dots ,x_m,y_1,y_2,\dots ,y_m\right]}^T, i=\mathrm{1,2,}\dots ,L$$ E1

Then the training image set $$\Omega =(X_1,X_2,\dots ,X_L)$$ .

The shape and pose (include position, scale and rotation) of the object in each training image is different, so the shapes should be aligned to filter out the pose effects. We will first explain how to align two shapes and then extend it to a set of shapes.

If X and X’ are two shape vectors, the goal is to minimize the square sum of corresponding landmarks after alignment using similarity transformation technique. That is to minimize E:

$$E={\left|T(X)-X^{\prime }\right|}^2$$ E2

where

$$T\left(\begin{array}{c}x\\ y\end{array}\right)=\left(\begin{array}{cc}u& -v\\ v& u\end{array}\right)\left(\begin{array}{c}x\\ y\end{array}\right)+\left(\begin{array}{c}{t}_x\\ t_y\end{array}\right)$$ E3

Choose three kinds of transformation factors: scale factor-s, rotation factor-θ and translation factor-t. And $$u=s\cdot \cos \theta$$ , $$v=s\cdot \sin \theta$$ . Hence $$s^2=u^2+v^2$$ $$\theta ={\tan }^{-1}(v/u)$$ . Then equation (2) can be rewritten as:

$$E(u,v,t_x,t_y)={\left|T(X)-X^{\prime }\right|}^2={\displaystyle \sum _{i=1}^m(u{x}_i-v{y}_i+t_x-{x^{\prime }}_i}{)}^2+{(v{x}_i+u{y}_i+t_y-{y^{\prime }}_i)}^2$$ E4

To solve the minimum value of equation (4) is equivalent to let the partial derivative equals to zero. Then the transformation parameters are:

$$u=\frac{(x\cdot x^{\prime })}{|x{|}^2},v=\frac{{\displaystyle \sum _{i=1}^m(x_i{y^{\prime }}_i-y_i{x^{\prime }}_i)}}{|x{|}^2},t_x=\frac{1}{m}{\displaystyle \sum _{i=1}^m{x^{\prime }}_i},t_y=\frac{1}{m}{\displaystyle \sum _{i=1}^m{y^{\prime }}_i}$$ E5

The alignment of a set of shapes can be processed by iterative approach suggested by [Bookstein, 1996]. The detailed process is shown in Table 1.

So the new training image set $$\widehat{\Omega }=({\widehat{X}}_1,{\widehat{X}}_2,\dots ,{\widehat{X}}_L)$$ . The mean shape is calculated by:

Step 1: Choose one shape as the initial estimate of the mean shape.

Step 2: Scale the mean shape so that$$\left|{\overline{X}}_0\right|=1$$.

Step 3: Align all the remaining shapes to the mean shape using the method described above.

Step 4: Re-estimate the mean shape from aligned shapes.

Step 5: Scale the new mean shape so that$$\left|{\overline{X}}_{new}\right|=1$$.

Step 6: If the new calculated mean shape doesn’t change significantly, convergence is declared; else return to step 3.

Table 1.

The process of aligning a set of shapes

$$\overline{X}=\frac{1}{L}{\displaystyle \sum _{i=1}^L{\widehat{X}}_i}$$ E6

The covariance matrix can thus be given as:

$${\displaystyle \sum {}_s=}\frac{1}{L-1}{\displaystyle \sum _{i=1}^L({\widehat{X}}_i-\overline{X}){({\widehat{X}}_i-\overline{X})}^T}$$ E7

Calculate the eigenvalue $${\lambda }_{s,i}$$ of $${\displaystyle \sum {}_s}$$ and the corresponding eigenvector $${\eta }_{s,i}$$

$${\displaystyle \sum {\eta }_{s,i}}={\lambda }_{s,i}{\eta }_{s,i}$$ E8

Sort all the eigenvalues in descending order:

$${\lambda }_{s,i}\ge {\lambda }_{s,i+1},i=\mathrm{1,2,}\dots \mathrm{,2}m-1$$ E9

The corresponding set of eigenvectors is $${{\rm H}}_s=[{\eta }_{s\mathrm{,1}},{\eta }_{s\mathrm{,2}},\dots ,{\eta }_{s\mathrm{,2}m}]$$

To reduce the computational complexity, the principal component analysis (PCA) method is used to reduce the space dimensions. Choose t largest eigenvalues which satisfy the following condition:

$${\displaystyle \sum _{i=1}^t{\lambda }_{s,i}}\ge 0.98\left({\displaystyle \sum _{i=1}^{2m}{\lambda }_{s,i}}\right)$$ E10

Then any shape instance can be generated by deforming the mean shape by a linear combination of eigenvectors:

$$\widehat{X}\approx \overline{X}+{\Phi }_s{b}_s$$ E11

where $${\Phi }_s$$ is a matrix of dimension $$2m\times t$$ and contains $$t$$ eigenvectors corresponding to the largest eigenvalues, $${\Phi }_s=({\phi }_{s\mathrm{,1}}|{\phi }_{s\mathrm{,2}}|\dots |{\phi }_{s,t})$$ . $$b_s$$ is a $$t$$ dimensional vector and the eigenvectors are mutually orthogonal, so it can be given by:

$$b_s={\Phi }_{-1}^s(\widehat{X}-\overline{X})={\Phi }_s^T(\widehat{X}-\overline{X})$$ E12

By varying the elements of b_s , new shape instance can be generated using equation (11).

2.3. Establish texture model

To establish a complete appearance model, not only the shape but also the texture should be considered. We should firstly acquire the pixel information over the region covered by the landmarks on each training image. Then the image warping method is used to consistently collect the texture information between the landmarks on different training images. Finally establish the texture model using the PCA method.

The texture of a target object can be represented as:

$$g={[g_1,g_2,\dots ,g_n]}^T$$ E13

where n denotes the number of pixel samples over the object surface.

In an image, the target object may occupy a small region of the whole image. The pixel intensities and global changes in illumination are different in each training image. The most important information is the texture that can reflect the characteristic of the target object. Due to the number of pixels in different target region is different and it is difficult to acquire the accurate corresponding relationship between different images, the texture model can not be established directly. We need to obtain a texture vector with the same dimension and corresponding relationship. So we warp each example image so that its control points match the reference shape. The process of image warping is described as follows:

Firstly, apply Delaunay triangulation to the reference shape to obtain the reference mesh which is consisted of a set of triangles. We choose the mean shape as the reference shape.

Secondly, suppose $$v_1$$ , $$v_2$$ and $$v_3$$ are three vertices of a triangle in the example mesh, any internal point $$v$$ of the triangle can be written as:

$$v=v_1+\beta (v_2-v_1)+\gamma (v_3-v_1)=\alpha v_1+\beta v_2+\gamma v_3$$ E14

where $$\alpha +\beta +\gamma =1$$ . Constrain $$0\le \alpha ,\beta ,\gamma \le 1$$ because $$v$$ is inside the triangle. Given $$v={[x,y]}^T$$ , $$v_1={[x_1,y_1]}^T$$ $$v_2={[x_2,y_2]}^T$$ $$v_3={[x_3,y_3]}^T$$ , then $$\alpha ,\beta ,\gamma$$ can be calculated by:

$$\alpha =1-(\beta +\gamma )$$ E15

$$\beta =\frac{y{x}_3-{\text{x}}_{\text{1}}y-x_3{y}_1-y_{3}x+{\text{x}}_{\text{1}}{y}_3+\text{x}{y}_1}{{\text{-x}}_{\text{2}}{y}_3+{\text{x}}_{\text{2}}{y}_1+{\text{x}}_{\text{1}}{y}_3+{\text{x}}_{\text{3}}{y}_2-{\text{x}}_{\text{3}}{y}_1-{\text{x}}_{\text{1}}{y}_2}$$ E16

$$\gamma =\frac{\text{x}{y}_2-\text{x}{y}_1-x_1{y}_2-{\text{x}}_{\text{2}}y+{\text{x}}_{\text{2}}{y}_1+{\text{x}}_{\text{1}}y}{{\text{-x}}_{\text{2}}{y}_3+{\text{x}}_{\text{2}}{y}_1+{\text{x}}_{\text{1}}{y}_3+{\text{x}}_{\text{3}}{y}_2-{\text{x}}_{\text{3}}{y}_1-{\text{x}}_{\text{1}}{y}_2}$$ E17

Finally, the corresponding point v’ of the triangle in the reference mesh can be calculated by:

$$v^{\prime }={v^{\prime }}_1+\beta ({v^{\prime }}_2-{v^{\prime }}_1)+\gamma ({v^{\prime }}_3-{v^{\prime }}_1)=\alpha {v^{\prime }}_1+\beta {v^{\prime }}_2+\gamma {v^{\prime }}_3$$ E18

where $${v^{\prime }}_1$$ , $${v^{\prime }}_2$$ and $${v^{\prime }}_3$$ are the vertices of the corresponding triangle in the reference mesh.

After image warping process, each shape in the training set is warped to the reference shape and sampled. The influence from the global linear changes in pixel intensities is removed.

To reduce the effects that caused by global lighting variations, the example samples are normalized by applying a scaling a and offset b:

$$g^{\prime }=({g^{\prime }}_1,{g^{\prime }}_2{g^{\prime }}_3,\dots ,{g^{\prime }}_n)=(\frac{g_1-b}{a},\frac{g_2-b}{a},\frac{g_3-b}{a},\dots ,\frac{g_n-b}{a})$$ E19

where

$$b=\frac{1}{n}{\displaystyle \sum _{i=1}^n{g^{\prime }}_i},a=\sigma ,{\sigma }^2=\frac{1}{n}{{\displaystyle \sum _{i=1}^n({g^{\prime }}_i-b)}}^2$$ E20

The establishment of the texture model is identical to the establishment of the shape model which is also analyzed by PCA approach. The mean texture is calculated by:

$$\overline{g}=\frac{1}{L}{\displaystyle \sum _{i=1}^L{g^{\prime }}_i}$$ E21

The covariance matrix can thus be given as:

$${\displaystyle \sum {}_g=}\frac{1}{L-1}{\displaystyle \sum _{i=1}^L({g^{\prime }}_i-\overline{g}){({g^{\prime }}_i-\overline{g})}^T}$$ E22

Calculate the eigenvalue $${\lambda }_{g,i}$$ of $${\displaystyle \sum {}_g}$$ and the corresponding eigenvector $${\eta }_{g,i}$$ :

$${\displaystyle \sum {\eta }_{g,i}}={\lambda }_{g,i}{\eta }_{g,i}$$ E23

Sort all the eigenvalues in descending order:

$${\lambda }_{g,i}\ge {\lambda }_{g,i+1},i=\mathrm{1,2,}\dots \mathrm{,2}n-1$$ E24

The corresponding set of eigenvectors is $${{\rm H}}_g=[{\eta }_{g\mathrm{,1}},{\eta }_{g\mathrm{,2}},\dots ,{\eta }_{g\mathrm{,2}n}]$$

To reduce the computational complexity, the PCA method is used to reduce the space dimensions. Choose t largest eigenvalues which satisfy the following condition:

$${\displaystyle \sum _{i=1}^t{\lambda }_{g,i}}\ge 0.98\left({\displaystyle \sum _{i=1}^{2n}{\lambda }_{g,i}}\right)$$ E25

Then any texture instance can be generated by deforming the mean texture by a linear combination of eigenvectors:

$$g^{\prime }\approx \overline{g}+{\Phi }_g{b}_g$$ E26

where $${\Phi }_g$$ is a matrix of dimension $$2n\times t$$ and contains t eigenvectors corresponding to the largest eigenvalues, $${\Phi }_g=({\varphi }_{g\mathrm{,1}}|{\varphi }_{g\mathrm{,2}}|\dots |{\varphi }_{g,t})$$ . $$b_g$$ is a t dimensional vector and the eigenvectors are mutually orthogonal, so it can be given by:

$$b_g={\Phi }_g^{-1}(g^{\prime }-\overline{g})={\Phi }_g^T(g^{\prime }-\overline{g})$$ E27

By varying the elements of $$b_g$$ , new texture instance can be generated using equation (26).

2.4 Establish combined model

Since the shape model and texture model have been established, any input image can be represented using the shape parameter vector b_s and texture parameter vector b_g . Since there are some correlations between shape and texture variations, the new vector b can be generated by combining b_s and b_g :

$$b=\left(\begin{array}{c}{W}_s{b}_s\\ b_g\end{array}\right)=\left(\begin{array}{c}{W}_s{\Phi }_s^T(\widehat{X}-\overline{X})\\ {\Phi }_g^T(g^{\prime }-\overline{g})\end{array}\right)$$ E28

where W_s is a diagonal matrix which adjust the weighting between pixel distances and pixel intensities. W_s is calculated by:

$$W_s=rI=\left[\begin{array}{ccc}r& \cdots & 0\\ ⋮& \ddots & ⋮\\ 0& \cdots & r\end{array}\right]$$ E29

where r ² is the ratio of the total intensity variation to the total shape variation:

$$r=\frac{{\lambda }_g}{{\lambda }_s},{\lambda }_g={\displaystyle \sum {\lambda }_{g,i}},{\lambda }_s={\displaystyle \sum {\lambda }_{s,i}}$$ E30

Apply PCA on b, then

$$b={\Phi }_{c}c$$ E31

where $${\Phi }_c$$ is the eigenvectors of covariance matrix corresponding to b:

$${\Phi }_c=\left(\begin{array}{c}{\Phi }_{c,s}\\ {\Phi }_{c,g}\end{array}\right)$$ E32

c is a vector of appearance model parameters controlling both the shape and texture of the models.

Using the linear nature of the model, the combined model including shape X and texture g can be expressed as:

$$X=\overline{X}+{\Phi }_s{W}^{-1}{\Phi }_{c,s}c$$ E33

$$g=\overline{g}+{\Phi }_g{\Phi }_{c,g}c$$ E34

Then a new image can be synthesised using equation (33) and (34) for a given c.

2.5 Fitting AAM to the input image

Fitting a AAMs to an image is considered to be a problem of minimizing the error between the input image and closest model instance [Wang et al., 2007]:

$$\delta I=I_{image}-I_{\mathrm{mod}el}$$ E35

Figure 3.

Results of tracking hand using AAM fitting.

where $$I_{image}$$ is the texture vector of the input image, and $$I_{\mathrm{mod}el}$$ is the texture vector of the model instance. The goal is to adjust the appearance model parameters $$c$$ to minimize $$|\delta I{|}^2$$ . The simplest way is to construct a linear relationship:

$$\delta c=R\delta I$$ E36

where R can be computed by the following process:

• Suppose c ₀ is the model parameter of the current image, new model parameter c can be generated by variating c ₀:

• Generate new shape model X and normalized texture model g_m according to equation (33) and (34).

• Deform the current image and get the corresponding texture model g_i , the difference vector can be written as:

$$\delta g=g_i-g_m$$ E38

$${\delta }_g$$ will change along with the variation of $${\delta }_c$$ and the shape model $$X$$

The fitting procedure is shown in Table 2 and the object tracking results are shown in Figure 3 and Figure 4. From the tracking results we can see that when the camera moves around the scene, the tracking results are satisfying.

3. Factorization algorithm

The flowchart of the factorization algorithm is shown in Figure 5, and the detailed demonstration will be given in the following sections.

3.1. Basic knowledge

The 3D shape of the non-rigid object can be described as a key frame basis set $$S_1,S_2,\dots ,S_K$$ . Each key frame basis $$S_i$$ is a $$3\times P$$ matrix describing $$P$$ points. The 3D shape of a specific configuration is a linear combination of the basis set:

Figure 4.

Results of tracking face using AAM fitting.

Step 1: Generate the normalized texture vector$$g_m$$.

Step 2: Sample the image $$g_i$$ below the model shape.

Step 3: Evaluate the error vector$${\delta }_{g_0}=g_i-g_m$$.

Step 4: Evaluate the error$$E_0=\left|{\delta }_{g_0}\right|$$.

Step 5: Calculate the pose displacement$${\delta }_t=R_t{\delta }_{g_0}$$.

Step 6: Calculate the displacement in model parameters$${\delta }_c=R_c{\delta }_{g_0}$$.

Step 7: Set$$i=1$$.

Step 8: Update model parameters$$c=c-k_i{\delta }_c$$.

Step 9: Transform the shape to invert the $${\delta }_t$$ transformation.

Step 10: Repeat step 1-4 to form a new error$$E_i$$.

Step 11: If $$E_i{E}_0$$ set $$i=i+1$$ and go to step 8.

Step 12: Accept the new estimate.

where $$k={\left[\mathrm{1.5,0.5,0.125,0.0125,0.0625}\right]}^T$$ is the damping vector.

Table 2.

The procedure of fitting the input image to the model instance.

$$S={\displaystyle \sum _{i=1}^K{l}_i\cdot S_i} S,S_i\in R^{3\times P},l_i\in R$$ E39

where

$$S=\left[\begin{array}{cccc}{x}_1& x_2& \mathrm{...}& x_P\\ y_1& y_2& \mathrm{...}& y_P\\ z_1& z_2& \mathrm{...}& z_P\end{array}\right]$$ E40

.

Under a scaled orthographic projection, the P points of S are projected into 2D image points (u_i v_i ):

$$\left[\begin{array}{cccc}{u}_1& u_2& \mathrm{...}& u_P\\ v_1& v_2& \mathrm{...}& v_P\end{array}\right]=R\left({\displaystyle \sum _{i=1}^K{l}_i{S}_i}\right)+T$$ E41

Figure 5.

The flowchart of the factorization algorithm.

$$R=\left[\begin{array}{ccc}{r}_1& r_2& r_3\\ r_4& r_5& r_6\end{array}\right],T=\left[\begin{array}{ccc}{t}_1& t_2& t_3\end{array}\right]$$ E42

R contains the first two rows of the 3D camera rotation matrix. T is the camera translation matrix. As mentioned in [Tomasi & Kanade, 1992], we eliminate the camera translation matrix by subtracting the mean of all 2D points, and henceforth can assume that the 3D shape S is centred at the origin.

$$\left[\begin{array}{cccc}{u^{\prime }}_1& {u^{\prime }}_2& \mathrm{...}& {u^{\prime }}_P\\ {v^{\prime }}_1& {v^{\prime }}_2& \mathrm{...}& {v^{\prime }}_P\end{array}\right]=\left[\begin{array}{cccc}{u}_1-\overline{u}& u_2-\overline{u}& \mathrm{...}& u_P-\overline{u}\\ v_1-\overline{v}& v_2-\overline{v}& \mathrm{...}& v_P-\overline{v}\end{array}\right]$$ E43

where

$$\overline{u}={\displaystyle \sum _{i=1}^P{u}_i},\overline{v}={\displaystyle \sum _{i=1}^P{v}_i}$$ E44

Therefore, we can rewrite equation (40) as:

$$\left[\begin{array}{cccc}{u^{\prime }}_1& {u^{\prime }}_2& \mathrm{...}& {u^{\prime }}_P\\ {v^{\prime }}_1& {v^{\prime }}_2& \mathrm{...}& {v^{\prime }}_P\end{array}\right]=R\left({\displaystyle \sum _{i=1}^K{l}_i{S}_i}\right)$$ E45

Rewrite the linear combination in equation (43) as a matrix-matrix multiplication:

$$\left[\begin{array}{cccc}{u^{\prime }}_1& {u^{\prime }}_2& \mathrm{...}& {u^{\prime }}_P\\ {v^{\prime }}_1& {v^{\prime }}_2& \mathrm{...}& {v^{\prime }}_P\end{array}\right]=\left[\begin{array}{cccc}{l}_{1}R& l_{2}R& \mathrm{...}& l_{K}R\end{array}\right]\cdot \left[\begin{array}{c}{S}_1\\ S_2\\ \mathrm{...}\\ S_K\end{array}\right]$$ E46

The 2D image points in each frame can be obtained using AAM. The tracked 2D points in frame t can be denoted as $$\left(u_{ti},v_{ti}\right)$$ . The 2D tracking matrix of N frames can be written as:

$$W=\left[\begin{array}{cccc}{u^{\prime }}_{11}& {u^{\prime }}_{12}& \mathrm{...}& {u^{\prime }}_{1P}\\ {v^{\prime }}_{11}& {v^{\prime }}_{12}& \mathrm{...}& {v^{\prime }}_{1P}\\ {u^{\prime }}_{21}& {u^{\prime }}_{22}& \mathrm{...}& {u^{\prime }}_{2P}\\ {v^{\prime }}_{21}& {v^{\prime }}_{22}& \mathrm{...}& {v^{\prime }}_{2P}\\ \mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}\\ {u^{\prime }}_{N1}& {u^{\prime }}_{N2}& & {u^{\prime }}_{NP}\\ {v^{\prime }}_{N1}& {v^{\prime }}_{N2}& & {v^{\prime }}_{NP}\end{array}\right]$$ E47

Using equation (44) we can rewrite equation (45) as:

$$W=\left[\begin{array}{cccc}{l}_{11}{R}_1& l_{12}{R}_1& \cdots & l_{1K}{R}_1\\ l_{21}{R}_2& l_{22}{R}_2& \cdots & l_{2K}{R}_2\\ \cdots & \cdots & \cdots & \cdots \\ l_{N1}{R}_N& l_{N2}{R}_N& \cdots & l_{NK}{R}_N\end{array}\right]\cdot \left[\begin{array}{c}{S}_1\\ S_2\\ \cdots \\ S_K\end{array}\right]$$ E48

where R_t denotes the camera rotation of frame t. l_ti denotes the shape parameter l_i of frame t.

3.2 Solving configuration weights using factorization

Equation (46) shows that the 2D tracking matrix $$W$$ has rank $$3K$$ , and can be factored into the product of two matrixes: $$W=Q\cdot B$$ , where $$Q$$ is a $$2N\times 3K$$ matrix, $$B$$ is a $$3K\times P$$ matrix. $$Q$$ contains the camera rotation matrix $$R_t$$ and configuration weights $$l_{t1},l_{t2},\cdots ,l_{tK}$$ of each frame. $$B$$ contains the information of shape basis $$S_1,S_2,\cdots ,S_K$$ . The factorization can be done using singular value decomposition (SVD) method:

$$W^{2N\times P}=\tilde{U}\cdot \tilde{D}\cdot {\tilde{V}}^T={\tilde{Q}}^{2N\times 3K}\cdot {\tilde{B}}^{3K\times P}$$ E49

Then the camera rotation matrix R_t and shape basis weights l_ti of each frame can be extracted from the matrix $$\tilde{Q}$$

$$q_t=\left[l_{t1}{R}_t{l}_{t2}{R}_t\cdots l_{tK}{R}_t\right]=\left[\begin{array}{ccccccc}{l}_1{r}_1& l_1{r}_2& l_1{r}_3& \cdots & l_K{r}_1& l_K{r}_2& l_K{r}_3\\ l_1{r}_4& l_1{r}_5& l_1{r}_6& \cdots & l_K{r}_4& l_K{r}_5& l_K{r}_6\end{array}\right]$$ E50

Transform $$q_t$$ into a new matrix $${\tilde{q}}_t$$ :

$${\tilde{q}}_t=\left[\begin{array}{cccccc}{l}_1{r}_1& l_1{r}_2& l_1{r}_3& l_1{r}_4& l_1{r}_5& l_1{r}_6\\ l_2{r}_1& l_2{r}_2& l_2{r}_3& l_2{r}_4& l_2{r}_5& l_2{r}_6\\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ l_K{r}_1& l_K{r}_2& l_K{r}_3& l_K{r}_4& l_K{r}_5& l_K{r}_6\end{array}\right]=\left[\begin{array}{c}{l}_{t1}\\ l_{t2}\\ ⋮\\ l_{tK}\end{array}\right]\cdot [\begin{array}{cccccc}{r}_{t1}& r_{t2}& r_{t3}& r_{t4}& r_{t5}& r_{t6}\end{array}]$$ E51

here, the $${\tilde{q}}_t$$ can be factored using SVD method.

3.3 Solving true rotation matrix and shape basis

As mentioned in [Tomasi & Kanade, 1992], the matrix $${\tilde{R}}_t$$ is a linear transformation of the true rotation matrix $$R_t$$ . Likewise, $${\tilde{S}}_i$$ is a linear transformation of the true shape matrix $$S_i$$ :

$$R_t={\tilde{R}}_t\cdot G,S_i=G^{-1}\cdot {\tilde{S}}_i$$ E52

where $$G^{3\times 3}$$ is found by solving a nonlinear data-fitting problem. In each frame we need to constrain the rotation matrix to be orthonormal. The constraints of frame t are:

$$\left[\begin{array}{ccc}{r}_{t1}& r_{t2}& r_{t3}\end{array}\right]G{G}^{-1}{\left[\begin{array}{ccc}{r}_{t1}& r_{t2}& r_{t3}\end{array}\right]}^T=1$$ E53

$$\left[\begin{array}{ccc}{r}_{t4}& r_{t5}& r_{t6}\end{array}\right]G{G}^{-1}{\left[\begin{array}{ccc}{r}_{t4}& r_{t5}& r_{t6}\end{array}\right]}^T=1$$ E54

$$\left[\begin{array}{ccc}{r}_{t1}& r_{t2}& r_{t3}\end{array}\right]G{G}^{-1}{\left[\begin{array}{ccc}{r}_{t4}& r_{t5}& r_{t6}\end{array}\right]}^T=0$$ E55

In summary, given 2D tracking data $$W$$ , we can get the 3D shape basis $${\tilde{S}}_i$$ , camera rotation matrix $${\tilde{R}}_t$$ and configuration weights $$l_{ti}$$ of each training frame simultaneously using factorization method.

4. Non-rigid registration method

The 2D tracking data can be obtained using the AAM algorithm mentioned in section 2. With the use of the factorization method mentioned in section 3, we can acquire the 3D shape basis. The 3D shape is represented as 3D points in the world coordinates. Given the configuration weights, the 3D shape can be recovered by linear combination of the 3D shape basis. By projecting the 3D points to the 2D image with known camera rotation matrix (suppose the intrinsic camera matrix has been calculated), the estimated 2D points can be acquired. If the error between the 2D tracking data and the estimated 2D points is small enough, we can accept the configuration weights and the rotation matrix. Finally, the virtual object can be overlaid to the real scene using the camera rotation matrix.

The initial configuration weights and camera rotation matrix can not be precise. Optimization of the configuration weights should be done to minimize the error between the 2D tracking points detected by AAM and the estimated 2D points which is projected by the 3D shape points. This is a non-linear optimization problem which can be successfully solved by the optimization methods. Different with [Zhu et al., 2006], we use the Levenberg-Marquardt algorithm. Equation (54) shows the cost function.

$${{\displaystyle \sum _j\Vert s_j-{s^{\prime }}_j\Vert }}^2\to \min ={{\displaystyle \sum _{j=1}^N\Vert s_j-\left(R{\displaystyle \sum _{i=1}^K{l}_i{S}_i}\right)\Vert }}^2\to \min$$ E56

where $$s_j$$ is the 2D tracking data $$j$$ $${s^{\prime }}_j$$ is the projected point.

The procedure of non-rigid registration is shown in Table 3.

Step 1: Track the 2D points $$s_j$$ using AAM.

Step 2: Initialize the configuration weights$$l_i$$.

Step 3: Initialize the camera rotation matrix$$R$$.

Step 4: Calculate the 3D shape$$S={\displaystyle \sum _{i=1}^K{l}_i{S}_i}$$

Step 5: Project the 3D points $$S$$ to 2D image:$${s^{\prime }}_j=A\left[R|T\right]\cdot S$$

Step 6: Evaluate the projection error$$E={{\displaystyle \sum _{j=1}^N\Vert s_j-{s^{\prime }}_j\Vert }}^2$$

Step 7: If $$E$$ is not small enough, improve $$l_i$$ and$$R$$, then repeat step 4-6.

Step 8: Accept the configuration weights$$l_i$$, and the camera rotation matrix$$R$$

Step 9: Overlaid the virtual object to the scene.

Table 3.

The procedure of non-rigid registration

Figure 6.

Examples of augmented images using our method.

Furthermore, we should take the orthonormality of rotation matrix into consideration. The proposed method has been implemented in C using OpenGL and OpenCV on a DELL workstation (CPU 1.7G×2, RAM 1G).

In offline stage, we construct the AAM hand models using 7 training images which are manually labelled with 36 landmarks as shown in Figure 1. We establish the hand shape basis using 300 training frames which is captured with a CCD camera. In online stage, the 2D points are tracked using AAM algorithm, and then Levenberg-Marquardt algorithm is used to optimize the parameters. Our experiment results are shown in Figure 6. From the results we can see that the virtual teapot can be overlaid on the hand accurately when the camera moves around the scene.

5. Conclusion

In this chapter, a non-rigid registration method for augmented reality applications using AAM and factorization method is proposed. The process is divided into two stages: offline stage and online stage. In the offline stage, the 3D shape basis is constructed. To obtain the shape basis of the object, we firstly factorize the 2D data matrix tracked by the AAM into the product of two matrixes. One matrix contains the camera rotation matrix and the configuration weights, and the other matrix contains the shape basis. Then the rotation matrix and the configuration weights can be separated using SVD method. Finally the orthonormality of the rotation matrix should be the constraints to get the true rotation matrix and configuration weights. In online stage, the 3D pose parameters and the shape coefficients are estimated. The purpose is to minimize the error between the 2D tracking points detected by AAM and the estimated 2D points which is projected by the 3D shape points. The Levenberg-Marquardt method is used to solve this problem. The rotation matrix and the configuration weights are optimized. Some experiments have been conducted to validate that the proposed method is effective and useful for non-rigid registration in augmented reality applications.

Acknowledgments

This work is supported in part by National Natural Science Foundation of China with project No. 60903095; in part by Postdoctoral Science Foundation of China with project No. 20080440941.