# Deformable Density Matching for 3D Non-rigid Registration of Shapes

Deformable Density Matching for 3D Non-rigid Registration of Shapes

Arunabha S. Roy1 , Ajay Gopinath1 , and Anand Rangarajan2,

1

Imaging Technologies Laboratory, GE Global Research Center, Bangalore, India 2 Department of CISE, University of Florida, Gainesville, FL, USA

Abstract. There exists a large body of literature on shape matching and registration in medical image analysis. However, most of the previous work is focused on matching particular sets of features—point-sets, lines, curves and surfaces. In this work, we forsake speci?c geometric shape representations and instead seek probabilistic representations— speci?cally Gaussian mixture models—of shapes. We evaluate a closedform distance between two probabilistic shape representations for the general case where the mixture models di?er in variance and the number of components. We then cast non-rigid registration as a deformable density matching problem. In our approach, we take one mixture density onto another by deforming the component centroids via a thin-plate spline (TPS) and also minimizing the distance with respect to the variance parameters. We validate our approach on synthetic and 3D arterial tree data and evaluate it on 3D hippocampal shapes.

1

Introduction

The need for shape matching occurs in diverse sub-domains of medical image analysis. Whenever a biomedical image is segmented or parsed into a set of shapes, the need for shape analysis and comparison usually arises [1]. In brain mapping for example [2], we frequently require the comparison of cortical and subcortical structures such as the thalamus, putamen etc. extracted from subject neuroanatomical MRI images. Image databases often use shape features and here the need is to index and query the shape database. In MR angiography, the complex network of blood vessels in the brain can be represented as trees or graphs and need to be compared across subjects. And in cardiac applications, if heart chamber information is available and extracted as a set of shapes, the wall tracking problem requires us to solve for shape correspondences in the cardiac cycle [3]. The need for shape matching is followed by a need for good shape representations. When shape features are extracted from medical images, they can be represented using an entire gamut of representations—points, line segments, curves, surfaces, trees, graphs and hybrid representations. What should be noted here is that an inferential stage is present in any shape representation. That is,

A.R. was partially supported by NSF IIS 0307712 and NIH RO1 NS046812.

N. Ayache, S. Ourselin, A. Maeder (Eds.): MICCAI 2007, Part I, LNCS 4791, pp. 942–949, 2007. c Springer-Verlag Berlin Heidelberg 2007

Deformable Density Matching for 3D Non-rigid Registration of Shapes

943

the raw features extracted from the underlying images are then converted into one of the above mentioned shape representations. Consequently, once a shape representation is adopted, we are immediately faced with the problem of robustness when we seek to compare two shapes. Some of the raw features present in (or extracted from) one image may not be present in the other. A second problem we face is that there may be no or poor correspondences between the features. It is for this reason that most methods seek to ?t curves or surfaces to the data. Once such representations are ?tted, there is no need to seek for correspondences at the point feature level. However, methods that rely on ?tting curves and surfaces to the data face a more di?cult robustness problem. How do you match one shape consisting of 10 curves to another shape consisting of 15 curves? For these reasons, we elect to go the probabilistic route. Beginning with raw point features, we ?t probability density functions to the feature vectors [4]. Each feature set is converted into a probability density function. The advantage here is that we can now compare probability density functions rather than the original images or other feature representations extracted from the images. The robustness problem is alleviated since the density functions are compared at all locations in R3 and we are not faced with the problem of matching incommensurate entities such as one curve in one shape to two curves in the other. There is no correspondence problem since we are not conducting comparisons at the point feature level. Shape comparison by matching the density functions between shapes also has the advantage that the point feature counts (in the two feature sets) can di?er considerably. We summarize our new method as follows: i) We ?t Gaussian mixture models to point features extracted from the two images. A standard maximum likelihood expectation-maximization (EM) approach is used for this step. ii) We derive a closed-form distance between the two Gaussian mixture models by comparing the probability density functions at each point in R3 . iii) Since the problem of minimizing this distance w.r.t. non-rigid deformations is ill-posed, we add a thin-plate spline regularization term to the cost. iv) We use a standard conjugategradient (CG) optimization strategy to minimize the above objective function. It should be stressed that we minimize the objective function w.r.t. both the deformation parameters and the variance parameters of the Gaussian mixture model. The result is a deformable density matching (DDM) method which seeks to register two shapes by moving the density function parameters of one shape until the ?rst density closely approximates the other.

2

Previous Work

There exists an enormous literature on feature matching. We restrict our focus to methods that seek to ?t and/or match probabilistic representations of shape features. The joint clustering and matching (JCM) approach [5] begins as we do with ?tting Gaussian mixture models to feature point-sets. However, instead of minimizing a distance between the two density functions, they seek to augment

944

A.S. Roy, A. Gopinath, and A. Rangarajan

the mixture model objective by linking a di?eomorphism between the centroid parameters of the mixture model. Consequently, they are forced to keep two sets of cluster centers in correspondence which we do not require. The methods in [6] and in [7] seek to convert a point matching problem into an image matching problem which is somewhat similar to density matching. However, the methods make no attempt to ?t a probabilistic shape representation via maximum likelihood (or equivalent) to the features. Essentially, both methods convert the sparse feature set into dense images and then employ an image matching strategy. Furthermore, the method in [6] uses a deformation ?eld parametrization which has to be applied at each point in R3 and this is computationally expensive compared to our approach. The method in [7] does not use a deformation ?eld parametrization but the main di?erences are that they restrict their approach to point matching and make no attempt to ?t a density function to the data via maximum likelihood or minimize their distance w.r.t. the centroid and variance parameters. Instead, they apply a thin-plate spline (TPS) on the original, noisy data. Perhaps the method that is closest in spirit to our approach is the recent work in [4]. They minimize the Jensen-Shannon divergence between the feature point-sets w.r.t. a non-rigid deformation. The Jensen-Shannon divergence cannot be computed in closed form for a mixture model and is estimated from the data using a law of large numbers approach. In sharp contrast, our distance between the two densities is in closed form and we apply the deformation parametrization to the centroidal parameters of the Gaussian mixture model.

3

3.1

Theory

Maximum-Likelihood Model for Shape Representation

As mentioned in the Introduction, the ?rst step in our overall method is probabilistic shape representation based on the raw shape features. The notation used in this paper is as follows. The two sets of input shape (1) (2) features are denoted by {Xi ∈ Rd , i ∈ {1, . . . , N1 }} and {Xj ∈ Rd , j ∈ {1, . . . , N2 }}. The maximum likelihood approach assumes that the shape fea(1) (2) tures {Xi } and {Xj } are independent and identically distributed (i.i.d.). The features of shape X (1) are represented by a Gaussian mixture model [8]

K1

p(x|θ(1) ) =

a=1

(1) Ωa

1 ?1 exp{? (x ? μa )T Σa (x ? μa )} 2 (2π) |Σa |

d 2 1 2

1

(1)

(where x ∈ R3 ) and that of shape X (2) is represented by a second Gaussian (2) mixture model with parameter set θ(2) = {Ωα , να , Ξα }. Constraints on Ω (·) (·) (·) K1 are {Ωa > 0, a=1 Ωa = 1} where the superscript index can be either 1 or 2 corresponding to the two shapes respectively. Both probability density functions de?ne measures on location with x ∈ R3 . We see that the set of (1) model parameters for shape X (1) is θ(1) = {Ωa , μa , Σa , a ∈ {1, . . . , K1 }} and

Deformable Density Matching for 3D Non-rigid Registration of Shapes

(2) (1)

945

(2)

θ(2) = {Ωα , να , Ξα , α ∈ {1, . . . , K2 }} for shape X (2) . Since {Xi } and {Xj } are assumed to be i.i.d., the likelihood of the set of features of X (1) is 1 (1) (1) ?1 exp{? (Xi ? μa )T Σa (Xi ? μa )} 2 (2π) |Σa | i=1 a=1 (2) (1) where {Xi } is the set of instances of X (1) . For both shapes, we ?t the model parameters by minimizing the negative log-likelihood objective function of the above mixture model w.r.t. the model parameters. In the experiments, we spe(1) 1 cialize to the case where the occupancy probabilities are uniform (Ωa = K1 ) 2 2 and where we have one isotropic covariance matrix (Σ = σ I3 ), for the entire shape. This is done for reasons of computational e?ciency. The objective function for this reduced version is minimized over its parameter set ({μa }, σ) to obtain the model representation of the feature set. We use the well known expectation-maximization (EM) algorithm [8] for the above minimization. While ?tting mixture models is computationally di?cult in the general case, in our special case of uniform occupancy probabilities and isotropic covariances it is not as di?cult. Also, this computation is done o?line for all the shapes once we ?x the number of centroids (K1 and K2 ). Model selection for mixture models needs to be performed to ?x the number of centroids.

(1) p({Xi }|θ(1) ) N 1 K1

=

(1) Ωa

1

d 2

1 2

3.2

A Closed-Form Distance Measure Between Two Gaussian Mixtures

We now derive the distance measure between the two probabilistic shape representations. Since this distance measure can be derived in closed form for the most general Gaussian mixture model case, we present this below. (Other distance measures like Kullback-Leibler cannot be derived in closed form for the Gaussian mixture model.) We hasten to add that the general mixture model can be quite unwieldy in practice and that specializations of the sort considered in the paper—like isotropic covariances and uniform occupancy probabilities—are very useful. Below we derive the distance measure as the squared pointwise difference between the two Gaussian mixture models with parameter sets (θ(1) and θ(2) ) integrated over Rd (where d = 3). [We have dropped terms that do not depend on θ(2) since it is only θ(2) that is deformed during the optimization.] D[p(x|θ(1) ), p(x|θ(2) )] =

K1 K2

?

a=1 α=1

2 exp{? 2(σ21 2 ) ||μa ? να ||2 } +ξ K1 K2 (σ 2 + ξ 2 ) 2

3

Rd K2 K2

[p(x|θ(1) ) ? p(x|θ(2) )]2 dx ∝

1 exp{? 4ξ2 ||να ? νβ ||2 } 2 2 2 K2 ξ 3

3

+

α=1 β=1

. (3)

This is the ?nal expression for the distance function used in this paper. The number of centroids and the variances of the two models can be di?erent.

946

A.S. Roy, A. Gopinath, and A. Rangarajan

3.3

Deformable Density Matching

We now turn to the description of the deformation model. We assume that the parameters θ(1) of shape 1 are held ?xed and that the parameters θ(2) = ({να }, ξ) (comprising the centroids and variance) of shape 2 are deformed so that p(x|θ(2) ) approaches p(x|θ(1) ). We use the familiar thin-plate spline (TPS) deformation model [9] for the centroids. That is, the action of the deformation takes the centroids {να } to

2 new locations {?α = Aνa + β=1 K(vα , νβ )Q2 γβ } where A is the unknown ν (4 × 3) a?ne matrix, the TPS kernel K(να , νβ ) = ?||να ? νβ || in 3D, Q2 is the K2 × (K2 ? 4) part of the QR decomposition of ν (in homogeneous coordinates) and γ is the unknown (K2 ? 4) × 3 matrix of deformation parameters. To avoid degenerate solutions which include all permutations of ν when K1 = K2 , we add a deformation regularization term λ trace(γ T QT KQ2 γ) to the objective function 2 in (3) with λ being a regularization parameter. In addition a regularization term λA trace[(A ? I)T (A ? I)] is added to prevent re?ections and unphysical a?ne transformations. We use a standard nonlinear conjugate-gradient (CG) algorithm (with line search) on the a?ne and deformation parameters (A, γ) and a 1-D search on the variance ξ 2 . The variance parameter is updated separately from the remaining parameters and is not allowed to abruptly change.

def

K

4

4.1

Results

Synthetic Example: Sphere to Ellipsoid Density Matching

To showcase our approach of density matching, we begin with a synthetic example. We generate a sphere and an ellipsoid with 600 points each. We run a mixture model EM algorithm on the ellipsoid and sphere representing them with 120 and 60 centroids respectively. We then warp the ellipsoid using a Gaussian radial basis function (GRBF) using the 120 centroids as the centers for the GRBF. Subsequently, we run the deformable density matching (DDM) on the two synthetic datasets. The goal is to deform the sphere so that it matches the warped ellipsoid. Here the ?xed variance σ 2 was set to the mean of the cluster variances resulting from the EM algorithm (σ = 0.04) and a 1-D search was performed over the parameter ξ, obtaining ξ = 0.06 as the optimum. The initial overlay of the sphere and the warped ellipsoid, the overlay of the centroids after matching and the ?nal overlay of the warped sphere (using the recovered deformation parameters) are shown in Figure 1. Our results clearly demonstrate that the lack of correspondences at the point level as seen in the ?nal shape overlay in Figure 1 is no deterrent to recovering the shape of the deformed ellipsoid. 4.2 Validation of Recovered Deformation: Bending and Stretching

Our approach to validation is based on comparing the recovered deformation vectors against the true (synthetically generated) deformation vectors generated at

Deformable Density Matching for 3D Non-rigid Registration of Shapes

947

Fig. 1. Left: Initial overlay of sphere and warped ellipsoid, Middle: Cluster centers of the two shapes after density matching and Right: Overlay after registration

1

0.14 0.12 0.1 0.08 0.06 0.04 0.02 0 bend stretch 0 0.2 0.4 0.6 warp factors 0.8 1

0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0 0.2 0.4 0.6 08 0.2 0.4 0.8 0.6

Fig. 2. Left: Vessel tree structure; Right: Validation over 30 noise trials, 10 warps

a set of vessel tree structure feature points [Figure 2]. Since we are not performing point matching, we cannot validate our results by using prior correspondence information. We begin with a point-set X consisting of about 200 points. A deformation ?eld is applied to the point-set X to obtain a deformed point-set Y . We remove 50% of the points in Y to get a reduced set and this is done so that the two mixture densities have a large discrepancy in the number of centroids. Then we match X to the reduced Y using deformable density matching using 60 and 30 centroids for the two sets respectively. The recovered TPS parameters are used to compute the estimated displacement vector at each point in X. We compare the estimated displacement to the true displacement by separately plotting the stretching and bending components of the mean-squared displacement uT ui ?i 1 1 1 u u u error: N2 i ||?i ? ui ||2 = N2 (||?i || ? ||ui ||)2 + N2 i 2||?i ||||ui ||(1 ? ||?i ||||ui || ) u where ui is the estimated deformation, ui is the true deformation and we have ? separated the total error into stretching (?rst term) and bending (second term) components. We executed 30 noise trials at di?erent TPS warp factors w ranging from 0.1 to 1 in steps of 0.1. The TPS warp factor speci?es that the TPS coe?cients be uniformly generated from the interval [?w, w]. Higher the warp factor, greater the degree of deformation. The medians of the two error measures are shown in Figure 2, illustrating the di?erence between bending and stretching. The bending errors are more sensitive and increase more rapidly with w.

bending and stretching validation

0.9

948

A.S. Roy, A. Gopinath, and A. Rangarajan

(a) Initial overlay of LATL (b) Final overlay of LATL(c) Final overlay of LATL hippocampal datasets hippocampal clusters hippocampal datasets Fig. 3.

(a) Initial overlay of RATL (b) Final overlay of RATL hippocampal datasets hippocampal clusters

(c) Final overlay of RATL hippocampal datasets Fig. 4.

4.3

Evaluation on 3D Hippocampal Datasets

The general problem of building a hippocampal atlas is a clinical application requiring multiple registrations of shapes from di?erent patients. Here we use deformable density matching to register 3D hippocampal point-set pairs from a database of 60 cases. We showcase results on two pairs of patients scheduled for left and right anterior temporal lobectomy (LATL and RATL) respectively. The point-sets consisting of a few hundred points each are clustered into about 80 centroids, and we set λ = 0.01. After matching, we overlay the two point-sets by warping one set onto the other using the recovered TPS parameters. The initial overlay of the data, the ?nal overlay of the centroids and the ?nal overlay of the

Deformable Density Matching for 3D Non-rigid Registration of Shapes

949

warped datasets are shown in Figures 3 and 4 for LATL and RATL respectively. In both cases, the ?nal overlay shows that we have achieved good registration. We are currently looking at entropic registration measures for more objective evaluation of all pairwise registrations across our database.

5

Discussion

In summary: i) We used a maximum likelihood EM approach to ?t centroids and variances to raw feature data. ii) We determined a closed-form distance between two Gaussian mixture models. iii) We recovered a TPS deformation of the centroids of one shape while also allowing the shape variance parameter to change in order to best match the two shape representations. iv) Finally, we applied the recovered TPS deformation to the original data and showed that shape matching was possible even though there were few correspondences at the point feature level. Our goal was to obtain a robust shape distance with very few free parameters and with the deformation parameters only a?ecting the model parameters and not speci?ed over R3 . The TPS model only warps the centroids and the variance parameters are allowed to move until the best ?t is achieved. There are four free parameters in our model: K1 , K2 and the regularization parameters λ, λA . Since the process of ?tting probabilistic representations is o?ine, we can take recourse to model selection in order to ?t K1 and K2 . The regularization parameters can be learned via a Bayesian MAP procedure. There also do not appear to be any technical barriers to generalizing the TPS deformation parametrization to a di?eomorphism and optimizing over the occupancy probabilities as well.

References

1. Davies, R.H., Twining, C., Cootes, T.F., Taylor, C.J.: A minimum description length approach to statistical shape modelling. IEEE Trans. Med. Imag. 21, 525–537 (2002) 2. Toga, A., Mazziotta, J.: Brain Mapping: The Methods. Academic Press, London (1996) 3. Tagare, H.: Shape-based nonrigid correspondence with application to heart motion analysis. IEEE Transactions on Medical Imaging 18(7), 570–579 (1999) 4. Wang, F., Vemuri, B.C., Rangarajan, A., Schmalfuss, I.M., Eisenschenck, S.J.: Simultaneous registration of multiple point-sets and atlas construction. In: Leonardis, A., Bischof, H., Pinz, A. (eds.) ECCV 2006. LNCS, vol. 3953, pp. 551–563. Springer, Heidelberg (2006) 5. Guo, H., Rangarajan, A., Joshi, S.C.: 3-D di?eomorphic shape registration on hippocampal data sets. In: Duncan, J.S., Gerig, G. (eds.) MICCAI 2005. LNCS, vol. 3750, pp. 984–991. Springer, Heidelberg (2005) 6. Glaunes, J., Trouve, A., Younes, L.: Di?eomorphic matching of distributions: A new approach for unlabelled point-sets and sub-manifolds matching. In: IEEE Conf. on Computer Vision and Pattern Recognition, vol. 2, pp. 712–718 (2004) 7. Jian, B., Vemuri, B.C.: A robust algorithm for point set registration using mixture of Gaussians. In: Intl. Conf. on Computer Vision (ICCV), pp. 1246–1251 (2006) 8. McLachlan, G.J., Basford, K.E.: Mixture models: inference and applications to clustering. Marcel Dekker, New York (1988) 9. Wahba, G.: Spline models for observational data. SIAM, Philadelphia, PA (1990)