Active Appearance Models: Theory and Cases

During the six months master thesis period, a paper was prepared and submitted to the 9th Danish Conference on Pattern Recognition and Image Analysis (DANKOMB). As documentation of the workload herein, the paper is reprinted below in one-column format. Since the paper was prepared in the middle of the thesis term, results may not be directly comparable to those reported in this thesis, due to changes in training set, implementation, evaluation methods and such. The nomenclature deviates also slightly from this thesis.





Active Appearance Models: Theory and Cases
M. B. Stegmann ^1,20.1, R. Fisker 1, B. K. Ersbøll 1,
H. H. Thodberg 2, L. Hyldstrup 3


1Department of Mathematical Modelling
Technical University of Denmark
DTU Building 321, DK-2800 Lyngby, Denmark


2Pronosco A/S, Kohavevej 5, DK-2950 Vedbæk, Denmark
3H:S Hvidovre Hospital, Kettegård Allé 30, DK-2650 Hvidovre, Denmark

Abstract

In this paper, we present a general approach towards image segmentation using the deformable model Active Appearance Model (AAM) as proposed by Cootes et al. A priori knowledge is learned through observation of shape and texture variation in a training set and is used to obtain a compact object class description, which can be used to rapidly search images for new object instances. An overview of the theory behind AAMs is given followed by an improved initialization scheme, thus making the AAMs fully automated. Finally, two cases are presented. It is demonstrated that AAMs can successfully segment bone structures in radiographs of human hands and structures of the human heart in 2D extracts of 4D cardiovascular magnetic resonance images. The observed mean point location accuracy was 1.0 and 1.3 pixels, respectively.
Keywords: Deformable Models, Snakes, Principal Component Analysis, Shape Analysis, Non-Rigid Object Segmentation, Initialization, Metacarpal Radiographs, Cardiovascular Magnetic Resonance Imaging.

Introduction

In recent years, the model-based approach towards image interpretation named deformable models has proven very successful. This is especially true in the case of images containing objects with large variability. Among the earliest and most well known deformable models is the Active Contour Model - known as Snakes proposed by Kass et al. [46]. Snakes represent objects as a set of outline landmarks upon which a correlation structure is forced to constrain local shape changes. In order to improve specificity, many attempts at hand crafting a priori knowledge into a deformable model have been carried out. These include Yuille's et al. [73] parameterization of a human eye using ellipsis and arcs. In a more general approach, while preserving specificity Cootes et al. [15] proposed the Active Shape Models (ASM) where shape variability is learned through observation. In practice this is accomplished by a training set of annotated examples followed by a Procrustes analysis combined with a principal component analysis. A direct extension of the ASM approach has lead to the Active Appearance Models [10]. Besides shape information, the textual information, i.e. the pixel intensities across the object, is included into the model. AAMs has been further developed in [14,22,13]. Quite similar to AAMs and developed in parallel herewith, Sclaroff & Isidoro suggested the Active Blob approach [58,43]. Active Blobs is a real-time tracking technique, which captures shape and textual information from a prototype image using a finite element model (FEM) to model shape variation. Compared to AAMs, Active Blobs deform a static texture, whereas AAMs change both texture and shape during the optimization. For further information on deformable models, refer to the surveys given in [44,51].

Active Appearance Models

Below we describe the outline of the Active Appearance Model approach. AAMs distinguish themselves in the sense that segmentation can be carried out using the approach as a black box. We need only provide with domain knowledge in the form of a training set annotated by specialists (e.g. radiologists etc.). Described is the training of the model, the modelling of shape and texture variation and the optimization of the model. Finally, an improved method for automated initialization of AAMs is devised. For a commented pictorial elaboration on the sections below - including the alignment process - refer to appendix A.

Shape & Landmarks

The first matter to clarify is: What do we actually understand by the term shape? This paper will adapt the definition by D.G. Kendall [20]:

Definition 6   $\hspace*{-1.0ex}\mathrm{:}$    Shape is all the geometrical information that remains when location, scale and rotational effects are filtered out from an object.

The next question that naturally arises is: How should one describe a shape? In everyday conversation unknown shapes are often described by referring to well-known shapes - e.g. "Italy has the shape of a boot". Such descriptions can obviously not be utilized in an algorithmic framework. One way of representing shape is by locating a finite number of points on the outline. Consequently the concept of a landmark is adapted [20]:

Definition 7   $\hspace*{-1.0ex}\mathrm{:}$    A landmark is a point of correspondence on each object that matches between and within populations.

A mathematical representation of an n-point shape in k dimensions could be concatenating each dimension in a kn-vector. The vector representation for planar shapes would then be:

$$\mathbf{x} = ( x_1, x_2,\dots , x_n, y_1,y_2, \dots, y_n )^{\mathrm{T}}$$ (20.1)

Notice that the above representation does not contain any explicit information about the point connectivity.

Shape Formulation

A classical statistical method for dealing with redundancy in multivariate data - such as shapes - is the   linear orthogonal transformation;  principal component analysis (PCA). In our application for describing shape variation by PCA - a shape of n points is considered one data point in a 2n^th dimensional space. In practice the PCA is performed as an eigenanalysis of the  covariance matrix of the shapes aligned w.r.t. position, scale and rotation, i.e. the shape analysis is performed on the true shapes according to the definition. As shape metric in the alignment the Procrustes distance [35] is used. Other shape metrics such as the Hausdorff distance [42] could also be considered. Consequently it is assumed that the set of N shapes constitutes some ellipsoid structure of which the centroid - the mean shape - can be estimated as:

$$\mathbf{\overline{x}} = \frac{1}{N} \sum^N_{i=1} \mathbf{x}_i$$ (20.2)

The maximum likelihood estimate of the covariance matrix can thus be given as:

$$\mathbf{\Sigma} = \frac{1}{N} \sum^N_{i=1} (\mathbf{x}_i - \m... ...line{x}}) (\mathbf{x}_i - \mathbf{\overline{x}})^{\mathrm{T}}$$ (20.3)

The principal axis of the 2n^th dimensional shape ellipsoid are now given as the eigenvectors, $\mathbf{\Phi}_s$, of the covariance matrix.

$$\mathbf{\Sigma} \mathbf{\Phi}_s = \mathbf{\Phi}_s {\lambda\hspace*{-1.22ex}\lambda}_s$$ (20.4)

A new shape instance can then be generated by deforming the mean shape by a linear combination of eigenvectors, weighted by b_s, also called the modal deformation parameters.

$$\mathbf{x} = \mathbf{\overline{x}} + \mathbf{\Phi}_s \mathbf{b}_s$$ (20.5)

Essentially, the point or nodal representation of shape has now been transformed into a modal representation where modes are ordered according to their deformation energy - i.e. the percentage of variation that they explain. What remains is to determine how many modes to retain. This leads to a trade off between the accuracy and the compactness of the model. However, it is safe to consider small-scale variation as noise. It can be shown that the variance across the axis corresponding to the i^th eigenvalue equals the eigenvalue itself, $\lambda_i$. Thus to retain p percent of the variation in the training set, t modes can be chosen satisfying:

$$\sum_{i=1}^t \lambda_i \geq \frac{p}{100} \sum_{i=1}^{2n} \lambda_i$$ (20.6)

Texture Formulation

Contrary to the prevalent understanding of the term texture in the computer vision community, this concept will be used somewhat differently below. The main reason for this is that most literature on AAMs uses this definition of texture, probably due to the close resemblance of some of the AAM techniques to techniques in computer graphics. In computer graphics, the term texture relates directly to the pixels mapped upon virtual 2D and 3D surfaces. Thus, we derive the following definition:

  Definition 3: Texture is the pixel intensities across the object in question (if necessary after a suitable normalization).

A vector is chosen, as the mathematical representation of texture, where m denotes the number of pixel samples over the object surface:

$$\mathbf{g} = ( g_1, \dots, g_m)^{\mathrm{T}}$$ (20.7)

In the shape case, the data acquisition is straightforward because the landmarks in the shape vector constitute the data itself. In the texture case one needs a consistent method for collecting the texture information between the landmarks, i.e. an image warping function needs to be established. This can be done in several ways. Here, a piece-wise affine warp based on the Delaunay triangulation of the mean shape is used. Another, theoretically better, approach might be to use thin-plate splines as proposed by Bookstein [5]. For details on the Delaunay triangulation and image warping refer to [33,60]. Following the warp sampling of pixels, a photometric normalization of the g-vectors of the training set is done to avoid the influence from global linear changes in pixel intensities. Hereafter, the analysis is identical to that of the shapes. Hence a compact PCA representation is derived to deform the texture in a manner similar to what is observed in the training set:

$$\mathbf{g} = \mathbf{\overline{g}} + \mathbf{\Phi}_g \mathbf{b}_g$$ (20.8)

Where $\mathbf{\overline{g}}$ is the mean texture; $\mathbf{\Phi}_g$ represents the eigenvectors of the covariance matrix and finally b_g are the modal texture deformation parameters. Notice that there will always be far more dimensions in the samples than observations thus leading to rank deficiency in the covariance matrix. Hence, to efficiently compute the eigenvectors of the covariance matrix one must reduce the problem through use of the Eckart-Young theorem. Consult [14,65] or a textbook in statistics for the details.

Combined Model Formulation

To remove correlation between shape and texture model parameters - and to make the model representation more compact - a 3rd PCA is performed on the shape and texture PCA scores of the training set, b to obtain the combined model parameters, c:

b = Q c (20.9)

The PCA scores are easily obtained due to the linear nature of the model:

$$\mathbf{b} = \left( \begin{array}{c} \mathbf{W}_s \mathbf{b... ...T}} (\mathbf{g} - \mathbf{\overline{g}}) \end{array} \right)$$ (20.10)

- where a suitable weighting between pixel distances and pixel intensities is obtained through the diagonal matrix W_s. An alternative approach is to perform the two initial PCAs based on the correlation matrix as opposed to the covariance matrix. Now - using simple linear algebra - a complete model instance including shape, x and texture, g, is generated using the c-model parameters.

$$\mathbf{x} =\mathbf{\overline{x}} + \mathbf{\Phi}_s \mathbf{W}_s^{-1} \mathbf{Q}_s \mathbf{c}$$ (20.11)

$$\mathbf{g} =\mathbf{\overline{g}} + \mathbf{\Phi}_g \mathbf{Q}_g \mathbf{c}$$ (20.12)

Regarding the compression of the model parameters, one should notice that the rank of Q will never exceed the number of examples in the training set. Observe that another feasible method to obtain the combined model is to concatenate both shape points and texture information into one observation vector from the start and then perform PCA on the correlation matrix of these observations.

Optimization

In AAMs the search is treated as an optimization problem in which the difference between the synthesized object delivered by the AAM and an actual image is to be minimized. In this way by adjusting the AAM-parameters ( c and pose) the model can deform to fit the image in the best possible way. Though we have seen that the parameterization of the object class in question can be compacted markedly by the principal component analysis it is far from an easy task to optimize the system. This is not only computationally cumbersome but also theoretically challenging - optimization theory-wise - since it is not guaranteed that the search-hyperspace is smooth and convex. However, AAMs circumvent these potential problems in a rather untraditional fashion. The key observation is that each model search constitutes what we call a prototype search - the search path and the optimal model parameters are unique in each search where the final model configuration matches this configuration. These prototype searches can be performed at model building time; thus saving the computationally expensive high-dimensional optimization. Below is described how to collect these prototype searches and how to utilizes them into a run-time efficient model search of an image. It should be noticed that the Active Blobs approach is optimized using a method quite similar to that of AAMs named difference decomposition as introduced by Gleicher [34].

Solving Parameter Optimization Off-line

It is proposed that the spatial pattern in $\delta \mathbf{I}$ can predict the needed adjustments in the model and pose parameters to minimize the difference between the synthesized object delivered by the AAM and an actual image, $\delta \mathbf{I}$:

$$\delta \mathbf{I} = \mathbf{I}_{image} - \mathbf{I}_{model}$$ (20.13)

The simplest model we can arrive at constitutes a linear relationship:

$$\delta \mathbf{c} = \mathbf{R} \delta \mathbf{I}$$ (20.14)

To determine a suitable R in equation (B.14), a set of experiments are conducted, the results of which are fed into a multivariate linear regression using principal component regression due to the dimensionality of the texture vectors. Each experiment displaces the parameters in question by a known amount and measuring the difference between the model and the image-part covered by the model.

  

Figure: Displacement plot for a series of y-pose parameter displacements. Actual displacement versus model prediction. Error bars are 1 std.dev.

As evaluation of the assumption of a linear relationship between the model and pose parameters and the observed texture differences, figure B.1 shows the actual and the mean predicted displacement from a number of displacements. The error bars correspond to one standard deviation. Hence, the optimization is performed as a set of iterations, where the linear model, in each iteration, predicts a set of changes in the pose and model parameters leading to a better model to image fit. Convergence is declared when an error measure is below a suitable threshold. As error measure, we use the squared L₂ norm of the texture difference, $\vert\delta \mathbf{g}\vert^2$. To gain a higher degree of robustness, one might consider using the Mahalanobis distance or a robust norm such as the Lorentzian error norm [58]. Fitness functions allowing for global non-linear transformations such as the mutual information [68,70] measure might also be considered.

Initialization

The optimization scheme described above is inherently sensitive to a good initialization. To accommodate this, we devise the following search-based scheme thus making the use of AAMs fully automated. The technique is somewhat inspired by the work of Cootes et al. [22]. The fact that the AAMs are self-contained or generative is exploited in the initialization - i.e. they can fully synthesize (near) photo-realistic objects of the class that they represent with regard to shape and textural appearance. Hence, the model, without any additional data, is used to perform the initialization. The idea is to use the inherent properties of the AAM-optimization - i.e. convergence within some range from the optimum. This is utilized to narrow an exhaustive search from a dense to sparse population of the hyperspace spanned by pose- and c-parameters. In other words, normal AAM-optimizations are performed sparsely over the image using perturbations of the model parameters. This has proven to be both feasible and robust. A set of relevant search configuration ranges is established and the sampling within this is done as sparsely as possible. Consider the graph given in figure B.1, which demonstrates that it should be safe to sample the y-parameter with a frequency of at least 10 pixels. One could also claim that as long the prediction preserves the right sign it is only a matter of sufficient iteration. To achieve sensitivity to pixel outliers we use the variance of the square difference vector between the model and the image as error measure:

$$e_{fit} = V[ \delta \mathbf{g}^2]\$$ (20.15)

As in the optimization this could easily be improved by using more elaborate error measures. In pseudo-code, the initialization scheme for detecting one object per image is:

1.

Set $e_{min} = \infty$ and m to a suitable low number (we use m=3)

2.

Obtain application specific search ranges within each parameter (e.g. $-\sigma \leq c_1 \leq \sigma$ etc.)

3.

Populate the space spanned by the ranges - as sparsely as the linear regression allows - by a set of sampling vectors $\mathbf{V} = \{ \mathbf{v}_1, \dots ,\mathbf{v}_n \}$.

4.

For each vector in V

5.

Do AAM optimization (max m iterations)

6.

Calulate the fit, e_fit, as given by (B.15)

7.

If e_fit < e_min Then e_min = e_fit, v_fit = v_n

8.

End

The vector v_fit will now hold the initial configuration. Notice that the application specific search ranges in step 2 is merely a help to increase initialization speed and robustness than it is a requirement. If nothing is known beforehand, step 2 is eliminated and an exhaustive search is performed. This approach can be accelerated substantially by searching in a multi-resolution (pyramidal) representation of the image.

Implementation

All experiments, illustrations etc. have been made using the Active Appearance Models Application Programmers Interface (AAM-API) developed in the C++ language and based on the Windows NT platform. This API will be released under the open source initiative in Autumn 2000^20.2, which means that other freely can download, use and elaborate on the AAM-API. The foundation for the AAM-API is the Intel Math Kernel Library for fast MMX implementation of BLAS, MS VisionSDK for image handling, ImageMagick for image I/O and finally DIVA for image processing and matrix handling. AAM-API performance compares to that of Cootes et al. As an example, the MRI optimizations took each 200 ms on average on a PII 350 MHz. Much effort has been put into providing documentation and educational features such as movies of the modes of variation, model search etc. Further info on AAMs, the AAM-API and full source code documentation can be obtained at the AAM-site.^20.3

Experimental Results

Segmentation in medical images has always posed a difficult problem due to the special image modalities (CT, MRI, PET etc.) and the large biological variability. To assess the performance of AAMs in such environments, our implementation has been tested upon radiographs of human hands and cardiac MRIs.

Radiographs of Metacarpals

Segmentation in radiographs (x-rays) pose a difficult problem due to large shape variability and inherent ambiguity of radiographs. This forms a suitable challenge. Other attempts to perform segmentation in radiographs include the work of Efford [25], where ASMs and other methods were used. Twenty radiographs of the human hand were obtained and three metacarpal bones were annotated using 50 points on each. The annotation of metacarpals 2,3 and 4 were concatenated into a 150-point model. To incorporate a more substantial texture contrast into the model, additional 150 points were placed along the normal of each model point, thus arriving at a 300-point shape model. The texture model consisted of approx. 10.000 pixels. Using 16 model parameters, 95% of the variation in the training set were explained. Automated initialization of the model followed by optimization reached a mean location accuracy of 0.98 pixels (point to associated border [9,17]) when testing on four unseen images with ground truth annotations. The mean texture error was approx. 7 gray levels (input images were in the byte range). Examples of initialization and optimization are given in the figures B.2-B.5. Notice a fairly good fit even in the distal (upper) and proximal (lower) end of the metacarpals where radiographs are rather ambiguous. To asses the performance within points, the mean point to point distance is plotted in figure B.6. Not surprisingly, problems arise in the distal and proximal end of the metacarpals due the large shape variability and the ambiguous nature of radiographs in regions of overlap.

  

Figure: Model border after automated initialization (cropped).

  

Figure: Optimized model border.

  

Figure: AAM after automated initialization (cropped).

  

Figure: Optimized AAM (cropped).

  

Figure: Mean point to point deviation from the ground truth annotation of each metacarpal. Low location accuracy is observed at the distal and proximal ends.

Cardiac MRIs

Another application in medical imaging is locating structures of 4D (space, time) cardiovascular magnetic resonance images. Temporal registration relative to the heart cycle has been done using ECG-triggered image acquisition. The pixel depth was 8 bits. An AAM has been built upon only four - spatially and temporally corresponding - 2D slices of four different hearts. The endocardial and epicardial contour have been annotated using 30 points each. The resolution of the 2D slices was 256x256 pixels resulting in a texture model of approx. 2600 pixels. The combined PCA explains 84 % of the variation in the training set using two model parameters. The cardiac AAM was then used to search in a unseen image spatially and temporally similar to those in the training set. The described initialization technique reached the result seen in figures B.7 and B.9. A final optimization reached a mean point accuracy of 1.7 pixels (point to associated border). The result can be seen in figures B.8 and B.10. The mean texture error was approx. 11 gray levels. By incorporating a model neighborhood similar to the metacarpal model the mean point accuracy was increased to 1.3 pixels. Future work on the cardiac AAM will include models built upon extracts differing spatially and temporally, thus leading to a somewhat unified AAM of a larger subspace from the original 4D. For a commented pictorial of the cardiac AAM, refer to appendix A. AAM segmentation of 2D cardiac MRIs has previously been done by Mitchell et al. [52]. A total of 102 images were used for the training set reaching a mean point accuracy of approx. 1 pixel on the endocardial and epicardial contour. Annotated structures were the right ventricle and endocardial and epicardial contours. The model was initialized manually.

  

Figure: Model border after automated initialization.

  

Figure: Optimized model border.

  

Figure: AAM after automated initialization (cropped).

  

Figure: Optimized AAM (cropped).

  

Figure: Original image (cropped).

Discussion & Conclusions

In this paper we have presented the basic theory of AAMs and devised a method providing sufficient initialization of AAMs. The performance of AAMs has been assessed on two different image modalities - i.e. radiographs and magnetic resonance images reaching a mean point location accuracy of 1.0 and 1.3 pixels, respectively. In both cases the location accuracy was noteworthy increased by adding a suitable neighborhood to the outer contours of the model, thus enhancing textual contrast. In the MRI case, we have shown that even with a training set as small as four examples, very good segmentation results can be obtained. This leads towards the straightforward assumption that the less variation observed in the object class in question, the smaller the training set one can allow. However - more than four images would probably still be desirable. The two cases stress the fact that the AAM approach is a example of a general vision technique that capture domain knowledge through observation. Contrary to this the bulk part of model-based vision techniques has hand crafted a priori knowledge by design and programming. Furthermore, we have experienced the AAM approach to be data-driven (non-parametric), self-contained and fast. We also notice that AAMs extend to multivariate imaging and higher spatial dimensions - i.e. into 3D models etc. More information on AAMs, papers, presentations, movies of model searches etc. can be obtained at the AAM-site at http://www.imm.dtu.dk/$\sim$aam/. A more comprehensive treatment of the work of M. B. Stegmann can be found in [65].

Acknowledgements

Cardiac MRIs were provided and annotated by M.D. Jens Christian Nilsson and M.D. Bjørn A. Grønning, H:S Hvidovre Hospital. M.Sc.Eng. Torben Lund provided the annotation tool. Dr. Rasmus Larsen, M.D. Jens Christian Nilsson and M.D. Bjørn A. Grønning is acknowledged for their very useful comments on the manuscript. The work of M. B. Stegmann is in part supported by a grant from Pronosco A/S.

Illustrated Cardiac AAM

  

Figure: Point cloud of four unaligned heart chamber annotations.

  

Figure: Point cloud of four aligned heart chamber annotations with mean shape fully drawn.

  

Figure: Correlation matrix of the four annotations. Observe the obvious point correlations.

  

Figure: Delanay triangulation of the mean shape.

  

Figure: Point variation of the four annotations; radius = $\sigma _x+\sigma _y$. Notice the large point variation to the lower left.

  

Figure: The first eigenvector plotted as displacement vectors. Notice that the large point variation observed in figure B.16 is point variation along the contour, which only contributes to a less compact model contrary to explaining actual shape variation.

  

Figure: Mean shape and shape deformed by the first eigenvector. Notice that this emphasizes the point above; that a lot of the deformation energy does not contribute to any actual shape changes.