1

Automatic Construction of Correspondences forTubular Surfaces

Toon Huysmans, Jan Sijbers, and Brigitte Verdonk

Abstract— Statistical shape modeling is an established tech-nique and is used for a variety of tasks in medical image pro-cessing, such as image segmentation and analysis. A challengingtask in the construction of a shape model is establishing agood correspondence across the set of training shapes. Especiallyfor shapes of cylindrical topology, very little work has beendone. This paper describes an automatic method to obtain acorrespondence for a set of cylindrical shapes. The method startsfrom an initial correspondence which is provided by cylindricalparameterization. The quality, measured in terms of the descrip-tion length, of the obtained correspondence is then improvedby deforming the parameterizations using cylindrical b-splinedeformations and by optimization of the spatial alignment of theshapes. In order to allow efficient gradient guided optimization,an analytic expression is provided for the gradient of this qualitymeasure with respect to the parameters of the parameterizationdeformation and the spatial alignment. A comparison is madebetween models obtained from the correspondences before andafter the optimization. The results show that, in comparisonwith parameterization based correspondences, this new methodestablishes correspondences that generate models with signifi-cantly increased performance in terms of reconstruction error,generalization ability, and specificity.

Index Terms— point correspondence problem, statistical shapemodels, tubular structures, minimum description length, imagesegmentation, image shape analysis.

I. INTRODUCTION

SHAPE correspondences and the derived statistical shapemodels have a wide range of applications in medical im-

age computing [1], [2]. They have been used to analyze shapedifferences between different classes of objects, for example,the lateral ventricles of schizophrenics versus a healthy group[3]. They have also been employed to gain more knowledgeabout the anatomical variability of certain organs or bones asfor example the human ear canal [4]. Such knowledge can inturn be used to reconstruct malformed, missing, or fracturedbone structures [5]. A widespread application of statisticalshape models is their use as prior knowledge in automaticimage segmentation [6]–[8]. The probability density functionof the shape is estimated from a set of manual segmentations.This knowledge is then used to guide the segmentation processof an unseen instance and to restrict the segmentation resultto the class of plausible shapes.

The major hurdle in the construction of a statistical shapemodel is establishing a dense correspondence over the surfacesof a large set of training shapes. These correspondencesshould be of high quality, i.e. the correspondence shouldmatch anatomically equivalent points over the surfaces. Ifthis requirement is not met, artificial modes of variation areintroduced into the shape model and this has a negative

effect on the performance of the model when used for imagesegmentation or interpretation [9].

For 2D shapes, a correspondence between the boundarycurves of the individual shapes is often defined by manuallandmarking [10]. Although this approach is feasible, it turnsout to be a time consuming and error prone task. In principle,the approach can be extended to 3D but it becomes highlyimpractical due to the large amount of landmarks that need tobe located and the increased level of difficulty in pin-pointingthem. Several approaches have been proposed to automate thislabor intensive procedure.

A relatively simple but effective approach is to establishthe sought-after correspondence by means of surface parame-terization [11]–[17]. Thereby, a one-to-one map is constructedbetween each surface of the set and some common, predefinedand usually mathematically simple parameter domain, such asa planar disc or a sphere. For each surface, the obtained one-to-one map associates a 2D parameter coordinate with eachpoint of the 3D surface. The surface-to-surface correspondenceis then defined by the assigned parameter values i.e., pointsbetween the surfaces correspond when they share the sameparameters.

Although the parameterization approach produces validcorrespondences, there is still room for improvement. Theparameterization of each shape in the training set is doneindependently of the other shapes, thus correlations betweenthe shapes are not taken into account. When an approachtakes this extra information into account, it can obtain bet-ter correspondences. Surface parameterization can provide agood initial correspondence for such a method. Davies et al.developed such a correspondence method in [18]. They use thedescription length of the derived shape model as a measure ofcorrespondence quality and optimize the correspondences withrespect to this criterion. To date, their minimum descriptionlength (MDL) approach is considered the method of choicefor the construction of correspondences [9].

Most of the aforementioned techniques focus on sets ofsurfaces of spherical topology since these are prevalent inbiomedical image processing, e.g. kidneys, liver, and brainventricles. Nonetheless, the developed principles can be trans-lated to surfaces of other topologies. More specifically, thispaper deals with the translation of these principles to sets ofsurfaces of cylindrical topology, such as the trachea, cochlearchannels, aortic aneurysms and the rectum. For such surfaces,the cylinder is the natural choice for the parameter domain. Aninitial correspondence is obtained by specialized cylindricalparameterization techniques. Following this, the correspon-dence is improved by an optimization framework that adopts

2

the MDL criterion as a measure of correspondence quality.The remainder of this paper is organized as follows. Related

work is discussed in Section I-A. An overview of the corre-spondence method is given in Section II. Parameterization ofsurfaces of cylindrical topology is discussed in Section III. Themeasure of correspondence quality, namely description length,is treated in Section IV. Section V details the procedure of howan initial correspondence is obtained and Section VI elaborateson the correspondence optimization. Experimental results ofapplying the method to a number of phantoms and real datasets are presented and discussed in Section VII. Section VIIIconcludes the paper.

A. Related Work

Parameterization of surfaces with disc-like topology ontoa planar convex region has been addressed in [11]. A para-meterization technique to parameterize surfaces of sphericaltopology onto the sphere has been developed by Brechbuhleret al. [12] and is known as SPHARM. This was later usedto model the shape of brain structures for segmentation [6]and analysis [3]. Conformal parameterization techniques forsurfaces of spherical topology were introduced in [19]–[21]. Amore efficient alternative to SPHARM, utilizing a progressivesurface representation, was provided by Praun and Hoppe [13].Based on this spherical parameterization technique, Huys-mans et al. have developed a cylindrical parameterizationtechnique to parameterize tubular surfaces onto the cylinderin a progressive way [14]. Conformal parameterization ofsurfaces of cylindrical topology was addressed in [22], [23]and [24]. Algorithms for more complex topologies (genus-n)have also been proposed [15]–[17] but, from a correspondenceoptimization point of view, these algorithms are suboptimalsince they rely on a heuristic or manually defined surfacechartification. In [25] Ricci flow is used to construct a seamlessperiodic tiling of a genus-n surface in the plane, sphere, orhyperbolic space and it is known as the universal coveringspace. The obtained map is angle preserving and can be usedto construct harmonic maps between surfaces of the sametopology [32].

In [26], Kottchef and Taylor used the determinant of thelandmark covariance matrix as an optimization objective.Based on their ideas, Davies et al. developed a minimumdescription length (MDL) formulation for the assessmentof correspondence quality for 2D curves [18], which waslater simplified by Thodberg [27]. The method for build-ing MDL correspondences has been extended to surfacesof spherical topology [28], which was later improved byHeimann et al. [29] in terms of computational efficiency.Horkeaw and Yang applied the MDL principle to surfacesof disc-like topology [30] and an extension to more complextopologies was obtained by cutting surfaces into topologicaldiscs prior to optimization [31]. This of course constrainsthe optimization along the cuts. In the recent work of Liet al. [32], a globally optimal map, in terms of harmonicenergy, is obtained between two surfaces sharing arbitrarycomplex topology. An extension of their method to populationsof surfaces and to other correspondence metrics, e.g. MDL,

is, however, not addressed. Recently introduced point-basedcorrespondence techniques can also handle arbitrary topology.In [34], Ferrarini et al. use self organising maps to obtaina pairwise correspondence between each population memberand a template. In [33], Cates et al. use particle systems tooptimize geometry sampling and groupwise correspondencewith respect to an information theoretic measure comparableto MDL. Point-based correspondence techniques are promisingbut problems occuring with highly convoluted surfaces stillneed to be addressed.

Little work has been done in shape modeling for cylindricalsurfaces. In [35], de Bruijne et al. proposed an improvedmodeling scheme for tubular objects. They used a manuallydetermined correspondence. Some of the previously discussedtechniques could be employed to build a correspondence fora set of cylindrical shapes. For example, the spherical MDLframework can be applied when the holes of the cylindricalshapes are closed but this can result in an invalid correspon-dence at the boundaries and this approach will not performwell for elongated surfaces [36]. The correspondence methodsthat rely on chartification [16], [31] could also be appliedto cylindrical surfaces but obviously the performance of thechartification heuristic will influence the final quality. Themethod proposed in this paper does not suffer from thesedrawbacks since, here, the description length minimizationmethod is translated specifically to treat sets of surfaces ofcylindrical topology.

II. METHOD OUTLINE

The method, treated in this paper, constructs a dense surfacecorrespondence together with an optimal spatial alignment foran arbitrary population of cylindrical surfaces. It proceeds intwo steps. First, a correspondence is derived from the surfaceparameterizations by alignment of the surfaces and parame-terizations. The result is refered to as the rigid correspon-dence. Then, the rigid correspondence is improved by applyinglocal, non-rigid, deformations to the parameterizations whilekeeping the surfaces optimally aligned. The finally obtainedcorrespondence is refered to as the non-rigid correspondence.An overview of these two steps can be found in Figure 1.

The input to the method is a set of ns triangle surfaces{M1, . . . ,Mns} of cylindrical topology. Each surface Mi isdefined by a tuple (Vi, Ti), where Vi is the set of nVi vertices{vi1, . . . ,vinVi} with vij ∈ R3 and Ti is the set of nTi triangles{ti1, . . . , tinT }.

The construction of the rigid correspondence is covered bythe flow chart in Figure 1(a) and presented in detail in SectionV. The rigid correspondence, denoted as {x1, . . . ,xns}, isobtained by parameterization of the surfacesMi, followed byrigid alignment of the surfaces and their parameter spaces. Thespatial transformations for the alignment of the surfaces aredenoted τ ◦(•|Φτ◦

i ) and the parameterization transformationsfor the parameter space alignments are denoted ρ◦(•|Φρ◦

i ).In order to avoid convergence to local minima, the optimalparameters Φτ◦

i and Φρ◦

i of these transformations are obtainedby a number of consecutive optimizations with respect tothe MDL correspondence quality criterion, denoted µ. The

3

construction pipeline for the rigid correspondence comprisesfollowing steps:

1) A cylindrical parameterization x◦i is constructed for eachsurfaceMi. It constitutes a one-to-one map between thesurface Mi and the open-ended cylinder C2

h of heighth.

2) A smooth b-spline approximation x◦i is constructedfor each surface. The b-spline representation results insmooth optimization objectives and provides a multi-resolution representation of the surface.

3) An initial spatial alignment of the surfaces is obtainedby alignment of their principal axes.

4) While keeping the spatial alignment fixed, an alignmentof the parameterizations is determined by optimizationwith respect to MDL.

5) Both the spatial and parameter space alignments areimproved by simultaneaous optimization with respect toMDL.

In order to obtain the rigid correspondence, the optimalspatial transformations and parameterization transformationsare applied to the b-spline approximations of each of theparameterized surfaces:

xi = τ ◦(•|Φτ◦i ) ◦ x◦i ◦ ρ◦(•|Φ

ρ◦

i ). (1)

The construction of non-rigid correspondence is covered bythe flow chart in Figure 1(b) and is presented in detail inSection VI. The non-rigid correspondence is obtained as animprovement of the rigid correspondence by applying a non-rigid b-spline transformation to the parameter space of eachsurface and simultaneously optimizing their spatial alignment.Again, the MDL correspondence criterion was used as theoptimization objective. For reasons of efficiency and in orderto avoid local minima, the optimization is done successivelyat a number of resolution levels. With each resolution level L,the grid size mL

u(0) × mLu(1) for the b-spline approximation

xLi of xi and the grid size nLu(0) × nL

u(1) for the b-splineparameterization deformation ρL(•|ΦρL

i ) is doubled, allowingfor a more detailed correspondence improvement. Also, thenumber of landmarks nLp used to calculate the shape model isincreased and the convergence tolerance becomes more strict.The spatial transformation τ (•|ΦτL

i ) is rigid and has optimalparameters ΦτL

i . The optimization of a resolution level L isinitialized with the result from the previous resolution levelL− 1. The final correspondence is obtained from the optimaltransformation parameters ΦτL

i and ΦρL

i of the last resolutionlevel:

xLi = τ (•|ΦτL

i ) ◦ xLi ◦ ρL(•|ΦρL

i ). (2)

From the final correspondence {x1, . . . , xns}, a map fromsurfaceMi to surfaceMj can be obtained by composition ofthe inverse of parameterization xi with the parameterizationxj , i.e. q = xj ◦ x−1

i (p), where p ∈Mi and the correspond-ing point q ∈Mj .

III. SURFACE REPRESENTATION

A. ParameterizationStarting from a cylindrical surfaceM defined by its vertices

V and triangles T , a parameterized version x ofM is obtained

(a) (b)

Fig. 2. A qualitative comparison of the harmonic parameterization methodwith the non-linear parameterization method. The parameterization of eachsurface is visualized by the blue and red lines which correspond to the iso-u(0)

and iso-u(1) lines of the cylinder, respectively. In (a), two parameterizationsof a surface with a large variation in cross-sectional diameter are shown.It is clearly visible that the harmonic parameterization (right) has largearea distortions at the bulge compared to the non-linear method (left), sothe non-linear method is preferred for this kind of surfaces. In (b), twoparameterizations of a surface with an approximately constant cross sectiondiameter are shown. For this surface, both methods can be used since theyapproximately perform equally well.

by assigning a unique pair of cylindrical coordinates to eachpoint of the surfaceM. Usually, the parameter coordinates areonly defined explicitly at the vertices of M and the extensionover the triangles is implied by barycentric interpolation ofthe parameter coordinates at the vertices. To be more precise,for surfaces of cylindrical topology, the parameterization x isa homeomorphic function from C2

h to the surface M, i.e.:

x : [0, 2π]× [0, h]→M⊂ R3

: u→ x(u),

where C2h denotes the open-ended two-dimensional cylinder

of length h with unit radius, which is parameterized by anangular coordinate u(0) and an axial coordinate u(1), i.e.u = (u(0), u(1)). For x to be a homeomorphism, it must be abijective, continuous function, and have a continuous inverse.If the topology of M is consistent, such a homeomorphismcan always be obtained although a solution is not unique. Thiscan be seen from the fact that the composition x ◦ ρ of anyautomorphism ρ of the cylinder with the parameterization xof the surface M, again is a valid parameterization of thesame surface. The particular solution that a parameterizationtechnique will propose is usually the result of the minimizationof an energy functional. Different functionals result in differentparameterizations, the quality of which depends on the actualapplication.

In order to obtain good correspondences between surfaces,a parameterization technique should create similar maps forsimilar surfaces. In addition, it is also desirable that it retainsrelative areas and angles as much as possible (i.e. distortion).When the parameterizations systematically suffer from largearea distortions, undersampling of parts of the surface canoccur in the final correspondence. For a tubular surface withan approximately constant width, the use of an harmonic para-meterization technique is recommended, e.g. [22], [23]. Theharmonic cylindrical parameterization is uniquely defined asthe solution of a system of linear equations. As a concequence,it is computationally very efficient. However, it fails to keep

4

cylindricalsurfaces

M1, . . . ,Mns

initialcorrespondencex1, . . . , xns

cylinderheight h

generatecylindrical

parameterizationsx◦1, . . . , x◦ns

b-spline surfacesfit x◦1 , . . . , x◦ns

apply spatialprincipal axes

alignment

τ◦(•|Φτ◦i )

converged?

yes

no

apply rigidparameterization

transform

ρ◦(•|Φρ◦i

)

L-BFGS

new Φρ◦i to

decrease µ

evaluatecorrespondence

quality µ

number oflandmarks np

evaluatecorrespondence

quality µ

new

Φτ◦i and Φρ◦

ito decrease µ

L-BFGS

apply transforms

ρ◦(•|Φρ◦i

) and

τ◦(•|Φτ◦i )

converged?no yes

(a) Rigid correspondence

initialcorrespondencex1, . . . , xns

b-splineapproximation

grid sizemL

u(0) ×mL

u(1)

b-spline surfaces

fit xL1 , . . . , xL

ns

proceed to nextresolution level

L ← L + 1

deformationb-splinegrid size

nL

u(0) × nL

u(1)

apply transformsτ(•|Φτ

i ) and

ρL(•|ΦρL

i)

new

Φτi and Φ

ρL

ito decrease µ

number of

landmarks nLp

evaluatecorrespondence

quality µ

convergencetolerances L-BFGS

stop multi-resolution?

no

yesconverged?

yes nooptimizedcorrespondencesx1, . . . , xns

(b) Non-rigid correspondence

Fig. 1. (a) Flow chart visualization of the rigid correspondence construction. First, each surface is equipped with a cylindrical parameterization andapproximated with a b-spline surface. Then, the surfaces are brought into a reference coordinate system by alignment of their principal axes. This isfollowed by the alignment of the parameter spaces of the surfaces by optimization w.r.t. the correspondence quality i.e., model description length. Finally,the rigid correspondence is obtained after a second optimization, where both the spatial and parameterization alignment parameters are set free. (b) Flowchart visualization of the non-rigid correspondence improvement. The correspondence is improved by simultaneously optimizing the parameters of the spatialalignments and the b-spline parameterization deformations w.r.t. the correspondence quality. The optimization is performed at multiple resolution levelssequentially, starting at the lowest level with coarse deformations and gradually adding more detail with each new level. Finally, the non-rigid correspondenceis obtained.

5

area distortions within acceptable bounds when large variationin cross-sectional diameter is present. For these cases, theprogressive non-linear cylindrical parameterization techniquefrom [14] is a better alternative. It allows control over thetrade-off between angle and area distortions, at the cost of anincreased computation time. See Figure 2 for a comparison ofthe two methods.

B. B-Spline Representation

Since the surfaceM is a piecewise linear surface, the partialderivatives of the surface coordinates with respect to the pa-rameter coordinates ( ∂x

∂u(0) and ∂x∂u(1) ) are discontinuous at the

triangle boundaries, excluding boundaries between coplanartriangles. This is undesirable since these partial derivativesare utilized in the gradient guided optimization of the corre-spondences. Using cubic b-splines, approximation of x willresult in a surface x that is C2 continuous within the b-splinepatches and C1 continuous at the patch boundaries. By varyingthe number of control points used for the approximation, amulti-resolution representation of the surface can be obtained.Moreover, evaluation of the b-spline representation x is muchfaster than evaluation of the triangle based representation x.This is because point location in a regular grid is far moreefficient than point location in a triangulation. A good b-splineapproximation of x can be achieved with a number of controlpoints much lower than the original number of vertices that xcomprises, resulting also in better memory efficiency.

The approximation uses a 2D tensor product b-spline sur-face. The cubic b-spline kernel, denoted β, is defined as in[37]:

β(u) =

16 (3|u|3 − 6u2 + 4), |u| ∈ [0, 1[16 (2− |u|)3, |u| ∈ [1, 2[0, |u| ∈ [2,∞[

. (3)

The b-spline surface is defined by a uniform grid of knotsK = {kij} positioned on the cylinder C2

h and a correspondinggrid of control points P = {pij} in R3. An example knotgrid is shown in Figure 4(a). The b-spline surface then hasthe following form:

x(u|K,P ) =mu(0)+1∑i=−1

mu(1)∑

j=−1

β

(u− kij

∆

)pij , (4)

where mu(0) is the number of knots in the u(0)-directionand mu(1) is the number of knots in u(1)-direction that arewithin the range [0, 2π]× [0, h]. The 2D cubic b-spline kernelβ(u) is separable, i.e. β(u) = β(u(0))β(u(1)) and the gridspacing is denoted as ∆ = ( 2π

nu(0)

, hnu(1)−1 ). The division

in the argument of the b-spline kernel is executed element-wise. In order to obtain a closed and smooth surface at theparameter boundary u(0) = 2π, the control points satisfy thethe following conditions:

p−1,j = pmu(0)−1,j

pmu(0) ,j = p0,j

pmu(0)+1,j = p1,j

(5)

Fig. 3. Approximation of a parameterized surface using a b-spline surfacewith an increasing number of control points (indicated as mu(0) ×mu(1) ).

Least squares fitting is used in order to find a set ofmu(0)(mu(1) + 2) control points that provides a good ap-proximation to the surface x. For this purpose, a set of mp

uniformly distributed parameter locations is chosen on thecylinder:Ump = {u1, . . . ,ump} ⊂ C2

h. Using these parameterlocations, the approximation error can be determined as thesum of the squared distances between the points on theoriginal surface and the points on the approximating surface atcorresponding parameter locations. The set of control pointsP for which this error is minimal, is regarded as the optimalset in a least squares sense:

P = arg minP

mp∑i=1

||x(ui|K,P )− x(ui)||2 (6)

The minimum is found as the solution of the system BP = X ,where X is a mp×3 matrix having surface points {x(ui)} asits rows, P is a mu(0)(mu(1) +2)×3 matrix having the controlpoints {pij} as its rows, and B is a mp ×mu(0)(mu(1) + 2)matrix. Note that Bkl is the contribution of the l-th controlpoint Pl to the approximation of the k-th surface point Xk.The solution is obtained efficiently from the normal equationsBTBP = BTX [38]. In Figure 3, a surface is shown togetherwith four cylindrical b-spline approximations obtained usingan increasing number of control points.

IV. CORRESPONDENCE QUALITY MEASURE

In this work, the quality of the correspondence of a set ofsurfaces is measured by the description length of the shapemodel that was built from the correspondence. The actualcomputation of the description length can be divided intothree steps. First, each surface is represented with a set ofcorresponding landmarks. Then, the point distribution model(PDM) is calculated from these landmarks. This PDM consistsof a mean surface, a set of shape modes and, the varianceexpressed by each of these modes. Finally, the descriptionlength is calculated from the obtained shape mode variances.The following three sections expound these steps.

A. Statistical Modeling

In order to build a statistical shape model or a pointdistribution model [1] for a set of ns surfaces {Mi}, acorrespondence needs to be established and the surfaces haveto be aligned in a common reference coordinate system.Suppose that {xi} are the parameterizations that express thiscorrespondence in the common reference coordinate system.

6

Then, the goal of statistical shape modeling is to capture theshape present in the set of surfaces with a probability densityfunction. Here, a distribution is assumed that is symmetricabout its mean namely a multivariate Gaussian distribution andthe actual distribution parameters are obtained using principalcomponents analysis (PCA) [39].

In this work, the computation of the PCA is done by meansof singular value decomposition (SVD). This is more efficientthan the traditional method where an eigenvalue decomposi-tion of the large covariance matrix is used. The SVD methodalso allows the computation of the partial derivatives of theshape mode variances w.r.t. the landmark positions, whichis important for the gradient based optimization. A matrixrepresentation of the set of surfaces {xi} is obtained by sam-pling each surface at a set of uniformly distributed cylindricalparameter locations Unp = {u1, . . . ,unp}. For each surface,the coordinates of the np landmarks are concatenated and a3np row vector xi, representing the surface xi, is obtained:

xi = [xi(u1) . . . xi(unp)]. (7)

The landmark matrix X is then obtained from the ns shapevectors as X = [xT1 . . . xTns ], resulting in a matrix ofdimensions 3np × ns. The mean shape vector x is computedas

x =1ns

ns∑i=1

xi (8)

and the row centered landmark matrix is obtained by subtract-ing this mean shape from each column of X , i.e.

Xc = [xT1 − xT . . . xTns − xT ]. (9)

Now, let the SVD of the centered landmark matrix be definedas

1√ns − 1

Xc = PSQT , (10)

where P is a 3np×3np orthonormal matrix containing the leftsingular vectors pj as its columns, S is a 3np × ns diagonalmatrix where the diagonal elements are the singular values σjin descending order, and Q is an ns×ns orthonormal matrixwith the right singular vectors qj as its columns. The ns×nssurface covariance matrix D is defined as

D =1

ns − 1XTc Xc = QS2QT (11)

and its SVD can be calculated efficiently. The first m columnsof P , denoted as Pm, contain the m shape modes pj and theyare obtained from the SVD of D as

Pm =1√

ns − 1XcQS

−1m , (12)

where Sm is the matrix that contains the first m rows of S.The first m corresponding shape mode variances λj are alsoobtained from Eq. (11) as the first m squared non-zero singularvalues σ2

j .Using the obtained shape modes pj and the corresponding

mode variances λj , a new shape instance x can be obtained

by adding a linear combination of the principal shape modesto the mean surface:

x = x+ns−1∑j=1

pjbj , (13)

where bj is the contribution of the j-th principal shape mode tox. Equation (13) defines a shape space spanned by the shapeparameters bj and with the mean shape as the origin. Thebounds on the shape parameters of the shape space are usuallychosen as a small multiple of the standard deviation of thepoint cloud along that direction, i.e. −3

√λj ≤ bj ≤ +3

√λj .

In what follows, Λ will denote the function that maps aset of corresponded surfaces to the mode variances of theirderived shape model:

(λ1, . . . , λns−1) = Λ(x1, . . . ,xns |Unp). (14)

B. Description Length

In Davies et al. [18], a correspondence measure for curvesand surfaces is introduced which is regarded as the currentstandard for correspondence optimization. Their measure isadopted here but in a simplified form. The original measureis based on the minimum description length principle: thesampled surfaces are coded in a message where the encodingis determined by the PCA model built from the correspon-dence. The total message length of the encoded surfaces,together with the encoded model, determine the quality of themodel and therefore also the quality of the correspondence.In this way, a trade-off is made between model complexityand goodness-of-fit. Over the years, the MDL measure hasbeen tuned and in this work the simplified MDL measure,introduced by Thodberg [27], is used. It is a function of theshape mode variances λj and is defined as follows:

µ(λ1, . . . , λns−1) =∑λi≥λc

(1 + log

λiλc

)+∑λi<λc

λiλc, (15)

The free parameter λc is set to be the expected noise variancein the data. The variation captured by all modes with aneigenvalue (variance) below λc is thus considered noise. Ascan be seen from the first and the second term in Equation(15) respectively, the benefit of decreasing normal modes islogarithmic while for noise modes it is constant. Furthermore,the quality measure µ goes to zero when all eigenvalues go tozero, i.e. it favors compact models. Also, both µ and its partialderivatives ∂µ

∂λiare continuous. This is an attractive property

for optimization. Note that shorter description length, i.e. alower value of µ, indicates better quality of correspondence.

C. Gradient of Description Length

The L-BFGS minimizer used in this work requires not onlythe value of the objective µ but also the gradient ∇µ. In thissection, the gradient with respect to the landmark positionsxij is derived. In [40], Ericsson and Kalle explain how toobtain the partial derivatives of the description length µ w.r.t.centered the landmark positions xcij . Their derivation is basedon a result obtained by Papadopoulo and Lourakis in [41].

7

This result is a simple expression, in function of the elementspik of P and qkj of Q, for the derivative of the singular valuesσk of Xc w.r.t. the matrix values xcij , namely:

∂σk∂xcij

= pikqkj . (16)

The partial derivatives of the description length w.r.t. the non-centered landmarks xij can be obtained as:

∂µ

∂xij=∑λk≥λc

1λk

∂λk∂xij

+∑λk<λc

1λc

∂λk∂xij

. (17)

The derivatives of the shape mode variances λk can be refinedinto:

∂λk∂xij

= 2σkpikqkj −2nsσkpik

∑

qk, (18)

where Equation (16) was used together with the fact that

∂xcı∂xij

=

1− 1

nsif ı = i and = j

− 1ns

if ı = i and 6= j

0 if ı 6= i

. (19)

Hereby the gradient of the description length w.r.t. thelandmarks are obtained. The gradient with respect to thetransformation parameters Φ can be obtained by multiplyingthe landmark gradient with the Jacobian of the function thatmaps the transformation parameters to the landmark positions.

V. RIGID CORRESPONDENCE

Establishing an initial correspondence for a set of surfaces{M1, . . . ,Mns} starts by parameterizing each surface ontothe cylinder. The parameterized piece-wise linear surfaces{x◦1, . . . ,x◦ns} are then approximated using b-splines, result-ing in the smooth surfaces {x◦1, . . . , x◦ns}. Details on this canbe found in Section III. The approximation uses a grid ofm◦u(0) × m◦

u(1) control points. It was observed that accurateapproximations are achieved with a grid size of the order32× h

2π32 for all surfaces considered in this work. An initialcorrespondence, denoted {x◦1, . . . , x◦ns}, is obtained from theb-spline surfaces after applying a spatial alignment τ ◦(•|Φτ◦

i )and parameter space rotation ρ◦(•|Φρ

◦

i ) to each of the surfacesx◦i . The optimal transformation parameters Φτ◦

i and Φρ◦

i aredetermined by optimization w.r.t. the model description lengthµ. The objective µ contains multiple local minima and there-fore suitable initialization for the transformation parametersis required. First the spatial transformation parameters areinitialized by principal axes alignment of the surfaces. Thisis followed by initialization of the parameterization rotationparameters by running a description length optimization for theparameterization rotations while keeping the spatial transfor-mation parameters fixed. Finally, starting from these initial pa-rameters, the full optimization of the description length, whereboth spatial transformation parameters and parameterizationrotation parameters are set free, is executed. This results inthe desired rigid correspondence {x◦1, . . . , x◦ns}. An overviewof establishing a rigid correspondence is given in the flowchart of Figure 1(a).

A. Spatial Alignment Initialization

The spatial transformation τ◦(•|Φτ◦) : R3 → R3, used forthe alignment of the surfaces, is a 3D rigid transformation. It iscomposed of a 3D rotation around the surface center followedby a 3D translation. The rotation is parameterized by a unitquaternion q = (w, qx, qy, qz) ∈ S3, where S3 is the three-sphere. In order to be of unit length, the quaternion shouldadhere to the following constraint:√

qqT = 1. (20)

Quaternion parameterization for rotation does not suffer fromthe singularities encountered with Euler angles. The transla-tion is parameterized by a 3d vector t = (tx, ty, tz). Thetransformation τ ◦ is thus controlled by seven parameters:Φτ◦ = (q, t).

Letx◦i (•|Φτ◦

i ) = τ ◦(•|Φτ◦i ) ◦ x◦i (21)

be a shorthand notation for the surface obtained after applyingthe spatial transformation τ ◦ to the surface x◦i , for givenparameters Φτ◦

i . Then, the parameters {Φτ◦1 , . . . ,Φτ◦

ns} arechosen so that the surfaces {x◦1(•|Φτ◦

1 ), . . . , x◦ns(•|Φτ◦ns)}

have their principal axes aligned with the axes of the referencecoordinate system. The translation vector ti for surface x◦icenters the surface at the origin of the coordinate system,i.e ti = −vi = − 1

nVi

∑nVij=1 v

ij where vi is the average of

the vertices vij of surface x◦i . The rotation for the surfaceis obtained using singular value decomposition. Let Vi =[(vi1 − vi)T . . . (vinVi − v

i)T ]T be the matrix that has thecentered vertices of Mi as its rows and let the singularvalue decomposition of the coordinate covariance matrix bedefined as 1

nVi−1VTi Vi = UiS

2iU

Ti . Then, Ui is the rotation

matrix that aligns the principal axes of the surface with thereference coordinate axes. Note that the rotation matrix Uiis not uniquely defined since the singular vectors are definedup to their sign. Thus there are eight possible rotations, fromwhich four can be eliminated since they produce a mirroredsurface. From the four remaining rotations, the one is chosenthat best matches a reference surface Mr in terms of thefollowing error:

np∑j=1

[D(x◦i (uj |Φτ◦

i ),Mr)]2, (22)

where D(p,M) measures the distance from p to the closestpoint onM. The quaternion qi that represents the best rotationtogether with the translation ti form the initialization for thespatial transformation parameter set, i.e. Φτ◦

i = (qi, ti).

B. Parameterization Alignment Initialization

The parameterization transformation ρ◦(•|Φρ◦), used toalign the parameterizations, is a parameter space rotation andit is controlled by a single parameter Φρ

◦ ∈ [0, 2π], definingthe angle of rotation. The transformation has the followingform:

ρ◦(•|Φρ◦) : C2

h 7→ C2h (23)

: u 7→ (u(0) + Φρ◦, u(1)). (24)

8

Similar to Equation (21), a shorthand notation for the surfaceobtained after applying the spatial and parameterization trans-formation, for given parameters Φρ◦

i and Φτ◦i , is as follows:

x◦i (•|Φρ◦

i ,Φτ◦i ) = τ ◦(•|Φτ◦

i ) ◦ x◦i ◦ ρ◦(•|Φρ◦

i ). (25)

The initial values for the parameterization rotation parameters{Φρ

◦

1 , . . . ,Φρ◦ns} are then obtained by solving the following

optimization problem:

arg minΦρ◦i ,∀i

µ ◦ Λ(. . . , x◦i (•|Φρ◦

i ,Φτ◦i ), . . . |Un◦). (26)

where the Φτ◦i are the initial spatial alignment parameters

obtained in the previous section.

C. Full Alignment

Now, starting from the initial spatial transformation param-eters {Φτ◦

1 , . . . ,Φτ◦ns} and the initial parameterization rotation

parameters {Φρ◦

1 , . . . ,Φρ◦ns}, obtained in the previous two

sections, the full description length minimization with rigidtransformations can be solved:

arg minΦρ◦i ,Φτ

◦i ,∀i

µ ◦ Λ(. . . , x◦i (•|Φρ◦

i ,Φτ◦i ), . . . |Un◦)

+ατ

ns

∑j

ητ (Φτ◦j ). (27)

Here, the regularization term ητ (Φτ◦j ) penalizes parameter sets

with a quaternion that violates Equation (20). The penaltyfor a quaternion q is measured as (1 −

√qqT )2. The reg-

ularization term attains its minimum when all quaternionsare in S3 and smoothly penalizes any deviation from this.In all the experiments a regularization factor ατ = 106 wasused. The parameters that solve Equation (27) are denoted{Φτ◦

1 , Φρ◦

1 , . . . , Φτ◦ns , Φ

ρ◦ns} and they provide the final rigid

correspondence {x◦1(•|Φρ◦

1 , Φτ◦1 ), . . . , x◦ns(•|Φ

ρ◦ns , Φ

τ◦ns)}.

VI. NON-RIGID CORRESPONDENCE

In the previous section, an rigid correspondence{x1, . . . ,xns} was produced by applying the optimalrigid spatial transformations and rigid parameterizationtransformations to the parameterized surfaces. Animprovement over this rigid correspondence can be obtainedby allowing local deformations of the parameterizations.Such local deformations can be realized by a non-rigidparameterization transformation. Here, a cylindrical b-splineparameterization deformation ρL is used, where L indicatesthe level of resolution of the deformation. The spatialtransformation τ is the same as in Section V. The optimalnon-rigid parameterization deformations, together withthe spatial transformation parameters, are determined byoptimization. The optimization is done at increasing levelsof resolution, sequentially, to avoid convergence to a localoptimum. The flow chart in Figure 1(b) gives an overview ofthe multi-resolution correspondence optimization.

At every resolution level L, the optimal transformation pa-rameters for that level are determined in the following manner:

First, each surface xi, obtained from the rigid correspondenceprocedure, is approximated with a b-spline surface, denotedxLi . The approximation uses a mL

u(0) ×mLu(1) grid of control

points. Here, mLu(0) = 3·2L−1 and mL

u(1) =⌊h2πm

Lu(0)

⌋is used,

thus the resolution is approximately isotropic and doubledfrom one level to the next. Each surface xLi is transformedaccording to the spatial transformation parameters ΦτL

i andparameterization transformation parameters ΦτL

i as follows:

xLi (•|ΦρL

i ,ΦτL

i ) = τ (•|ΦτL

i ) ◦ xLi ◦ ρL(•|ΦρL

i ). (28)

The level-L b-spline parameterization transformation ρL isdefined by a grid of nL

u(0) × nLu(1) control points. Similar tothe approximation b-spline, the transformation b-spline hasan isotropic resolution that is doubled at each new level, i.e.nLu(0) = 3 · 2L−1 and nL

u(1) =⌊h2πn

Lu(0)

⌋. Using the notation

from Equation (28), the optimal level-L transformation pa-rameters are determined by solving the following optimizationproblem:

arg minΦρ

L

i ,ΦτLi ,∀i

µ ◦ Λ(. . . , xLi (•|ΦρL

i ,ΦτL

i ), . . . |UnLp )

+ατ

ns

∑j

ητ (ΦτL

j ) +αρ

ns

∑j

ηρL

(ΦρL

j ). (29)

Here, ητ is the regularization for the spatial transformations,as defined in Section V-C, and ηρ

L

is the regularization forthe parameterization deformations in order to avoid overfitting.The experimentally determined regularization constants areατ = 106 and αρ = 0.2. The set of parameter coordinatesUnLp , used to estimate the shape covariance matrix, containsnLp = 250 · 4L−1 parameter locations. Thus, the number ofsamples per area increases fourfold with every new level.This mirrors the doubling of the resolution of the surfaceapproximation and parameterization transformation. Note thatthe optimal transformation parameters for the optimizationproblem of level L − 1 are used to initialize the parametersof current the level L. To initialize the b-spline parameteriza-tion transformation parameters ΦρL

j from ΦρL−1

j the b-splineupsampling technique from [42] is used. The initialization ofthe spatial transformations is trivial, i.e ΦτL

j = ΦτL−1

j . Inthis work, three levels of resolution were used and thus theoptimal transformation parameters {Φτ3

1 , Φρ3

1 , . . . , Φτ3

ns , Φρ3

ns}of the third level optimization problem are the final trans-formation parameters. These provide the final correspondence{x1, . . . , xns}.

In the following two sections, the actual form of theparameterization transformation will be detailed and a suitableregularizer is introduced.

A. Reparameterization Transformation

The parameterization space C2h is deformed using a para-

meterization transformation ρL(u|ΦρL), where L denotes thelevel of resolution. The transformation is an automorphismof the parameter space, i.e. ρL constitutes a continuous one-to-one map of C2

h. The space of possible reparameterizationsis spanned by the transformation parameters ΦρL . Differentbases can be used to represent a reparameterization function

9

[18], [28], [29], [43]. In this work, ρL is a 2d cubic b-spline function with knot positions on a regular grid. Sucha representation has a number of convenient properties: (1)cubic b-spline deformations are C2 continuous with respectto their parameters within the patches and C1 continuous atthe patch boundaries. This is required for efficient, gradientguided, optimization. (2) it has compact support which makesit fast to evaluate and allows local control. And, (3) it canbe used in a multi-resolution method by refining the grid thatcontrols the shape of the deformation.

The b-spline deformation function ρL is defined by a set ofknots KL = {κLij} and a set of control point displacementsΦρL = {δLij}:

ρL(u|ΦρL) =nLu(0)+1∑i=−1

nLu(1)∑

j=−1

β

(u− κLij

∆L

)(κLij + δLij) mod(0) 2π (30)

Where the modulo operator, mod(0), acts on the first coordi-nate of the parameter space and keeps the deformed parameterwithin the bounds [0, 2π]. β is the 2D separable cubic b-splinekernel from Equation (3). The knots and the correspondingcontrol points are arranged on a regular nL

u(0) × nLu(1) grid onthe cylinder C2

h. The spacing between the knots is denoted as∆L = ( 2π

nLu(0)

, hnLu(1)−1

).

The following constraints on the control point displacementsδij make sure that ρL is periodic and continuous at theparameter boundary u(0) = 2π:

δ−1,j = δnu(0)−1,j

δ0,j = δnu(0) ,j

δ1,j = δnu(0)+1,j

,∀j. (31)

In order to make ρL one-to-one along the boundaries of C2h,

the following constraints are also enforced:δ

(1)i,−1 = −δ(1)

i,1

δ(1)i,n

u(0)= −δ(1)

i,nu(0)−2

δ(1)i,0 = 0δ

(1)i,n

u(0)−1 = 0

,∀i. (32)

In this work, the same arrangement of knots KL is used forthe reparameterization of each surface xk. Only the controlpoint displacements ΦρL

k = {δLij,k} differ from surface tosurface. Figure 4(a) shows a 3 × 4 cylindrical grid of knotsin an unfolded view and Figure 4(b) shows a grid of controlpoints for that knot grid, where the control points are obtainedby adding the control point displacement to the knot location,i.e. κLij+δ

Lij . Observe that both the knots and the control points

reside on the parameter domain C2h. The b-spline function ρL

thus defines a reparameterization function of the cylindricalparameter domain C2

h.

B. Reparameterization Regularization

From Equation (30), it can be seen that the reparameteriza-tion transformation ρL is a linear combination of translated

versions of the cubic b-spline kernel β. As a result, thereparameterization transformation is C2 continuous within thepatches and C1 continuous at the patch boundaries. The in-herent smoothness of the b-spline transformation is convenientbut it does not avoid overfitting by the b-spline transform.Overfitting is perceived as an irregular local deformation andoccurs in regions where µ is insensitive to local deformations,e.g. when a b-spline kernel is not supported by any landmarksamples in the calculation of Λ.

To counter these irregularities, a regularization term is intro-duced, denoted ηρ. For b-splines, a number of regularizationterms have been used in the past [44], [45]. Here, a simpleregularizer is chosen that measures the Dirichlet energy ofthe parameterization displacement function. This displacementfunction, denoted % can be derived easily from Equation 30:

%(u|Φρ) =nu(0)+1∑i=−1

nu(1)∑

j=−1

β

(u− κij

∆

)δij (33)

The regularization term ηρ for the deformation is then definedas:

η(Φρ) =∫∫C2h

(∣∣∣∣∂%(u|Φρ)∂u(0)

∣∣∣∣2 +∣∣∣∣∂%(u|Φρ)

∂u(1)

∣∣∣∣2)

du (34)

In what follows, βiu(c) = β

(u(c)

∆(c) − i)

is used as a shorthandnotation and the prime mark symbol denotes the derivative.Substituting the b-spline transformation into Equation (34)results in:

η(Φρ) =∑ij

∑kl

δijδTkl(

1

∆(0)2

∫ 2π

0

βi′

u(0)βk′

u(0)du(0)

∫ h

0

βju(1)β

lu(1)du(1)

+1

∆(1)2

∫ 2π

0

βiu(0)βku(0)du(0)

∫ h

0

βj′

u(1)βl′

u(1)du(1)

)(35)

From this, the derivatives of the smoothness energy w.r.t. thecontrol point displacements are easily derived:

∂η

∂δ(0)kl

= 2∑ij

δ(0)ij(

1

∆(0)2

∫ 2π

0

βi′

u(0)βk′

u(0)du(0)

∫ h

0

βju(1)β

lu(1)du(1)

+1

∆(1)2

∫ 2π

0

βiu(0)βku(0)du(0)

∫ h

0

βj′

u(1)βl′

u(1)du(1)

)(36)

Equations (35) and (36) can be evaluated analytically and, asa result, they can be used efficiently with a gradient basedoptimization method.

VII. RESULTS AND DISCUSSION

A. Data Sets

In this work, six populations of surfaces were used totest the proposed correspondence method. Three of theseare phantom populations, each consisting of 30 surfaces: a

10

u(0)

u(1)

κ−1,−1

κ−1,0

κ−1,1

κ−1,2

κ−1,3

κ−1,4

κ0,−1

κ0,0

κ0,1

κ0,2

κ0,3

κ0,4

κ1,−1

κ1,0

κ1,1

κ1,2

κ1,3

κ1,4

κ2,−1

κ2,0

κ2,1

κ2,2

κ2,3

κ2,4

κ3,−1

κ3,0

κ3,1

κ3,2

κ3,3

κ3,4

κ4,−1

κ4,0

κ4,1

κ4,2

κ4,3

κ4,4

2π

h

(a) knot grid

u(0)

u(1)

κij

δij

2π

h

(b) control points

Fig. 4. (a) A 3×4 grid of knots on C2h. The knots κi,−1 and κi,4 ensure that the cubic b-spline can be evaluated over the whole domain of C2h (indicated asthe shaded area) and the knots κ−1,j , κ3,j , and κ4,j make the b-spline periodic in the u(0)-direction. (b) An example grid of b-spline control points overlayedon the knot grid. Together these define a cylindrical parameterization transformation. The shape of the transformation is controlled by the displacements δijfrom the corresponding knots κij .

(a) synthetic (b) real

Fig. 5. From each of the six considered surface populations one randomsurface is shown here. Each surface is textured with the iso-contours of itsparameter coordinates obtained after parameterization. The three populationsof synthetic surfaces, in clockwise order in (a) are cylinders, disks, and beams.The three populations of real surfaces, in clockwise order in (b) are clavicles,tracheas, and thrombi.

population of disks where the position of the disk on thecylinder is variable, a population of beams with varying widthand depth, and a population of cylinder-like bent surfaceswith elliptic cross-section where the width and height of thecross-section is variable together with the amount of bending.A sample surface of each of these populations is shown inFigure 5(a). There are also three populations of real, ct-scanned, surfaces: a population of 25 clavicles, a populationof 23 tracheas and a population of 50 aortic sections with athrombus. A sample surface of each of the real populations canbe found in Figure 5(b). All surfaces of these populations havecylindrical topology except the clavicles which are of sphericaltopology. The clavicles are made cylindrical by puncturingboth ends of the clavicle.

B. Performance Measures

In what follows, different correspondences will be con-structed for the above-mentioned populations. The quality ofthe established correspondences is evaluated by deriving aPCA-model from the correspondence and reporting perfor-mance measures for the obtained model. The performance ofa model is measured here by the compactness, reconstructionability, generalization ability and specificity of the model. Theperformance measures are reported for the full k-mode modeland all restricted m-mode, m < k, versions of the PCA-model. In a comparison, the correspondence having the bestperformance measures, for its derived model, is considered thebest correspondence.

Now, given a set of surfaces {M1, . . . ,Mns}, let thecorrespondence be denoted as {x1, . . . ,xns} and the derivedshape model as xm, where m is the number of modes of theshape model. Then the compactness of a model is measuredas the cumulative variance:

C(m) =m∑i=1

λi, (37)

where λi is the variance of the i-th shape mode. The re-construction ability indicates how good a model is able toreconstruct the surfaces that were used to build the model. Itis measured as the average approximation error after fittingthe model to each of the surfaces {M1, . . . ,Mns}:

R(m) =1ns

ns∑i=1

minΦτ ,b

D(τ (Φτ ) ◦ xm(b),Mi), (38)

where Φτ are the parameters of the rigid transformation τ ,b are model parameters, and D(x,y) measures the averageclosest point distance from surface y to surface x. The optimalparameters Φτ and b, resulting in the best model-to-surface

11

i= 1

i= 2

i= 3

x x+3biqλix − 3bi

qλi x x+3bi

qλix − 3bi

qλi

trachea clavicle thrombus

x x+3biqλix − 3bi

qλi

rigid

bspline

i = 1

i = 1

i = 2

i = 3

i = 4

i = 5

disk

x x+3biqλix − 3bi

qλi

Fig. 6. A visualization of the modes of models obtained with the approach of this paper. The color coding represents the frobenius norm of the landmarkcovariance matrix and is a measure of the local variability of the surface. On the left, the first three modes are shown for the three CT-scanned populations:trachea, clavicle, and thrombus. For each mode i, the average shape x is shown together with the positive and negative offset in the direction of mode i:x± 3bi

√λi. On the right, the modes of model for the disk population is shown. On top, the b-spline optimized model is shown. At the bottom, the model

of the rigid correspondence is shown. Clearly, the b-spline optimization dramatically improved the correspondence.

fit, are determined iteratively by alternately estimating Φτ andb in a least squares sense. The generalization ability of amodel determines how well the model generalizes to unseeninstances of the modeled class. It is measured as the averageapproximation error after fitting leave-one-out versions of themodel to the left out surfaces:

G(m) =1ns

ns∑i=1

minΦτ ,b

D(τ (Φτ ) ◦ xmi (b),Mi), (39)

where xmi is the m-mode model where the i-th surfacewas left out, i.e. it is built from the corresponded surfaces{x1, . . . ,xi−1,xi+1, . . . ,xns}. The model specificity mea-sures how much random samples, generated by the model,resemble the original surfaces:

S(m) =1nt

nt∑i=1

minj,Φτ

D(τ (Φτ ) ◦ xm(bmi ),Mj), (40)

where the bmi are random Gaussian model parameters forthe sample of the i-th trial and nt is the number of randomsamples used to estimate the specificity. Note that the modelperformance measures from [18] result in a bias towards theMDL-optimized models. The model performance measures inEq. (38), (39), and (40) do not suffer from this drawback.The compactness measure from Eq. (37) is closely relatedto the MDL-measure and therefore biased. It is, however,reported here because it containes important information ofhow a model captures the variation of a population.

MDLpopulation ns h ideal rigid non-rigid time (min.)

disks 30 2π 6.2 17.2 4.6 66beams 30 2π 9.6 9.8 9.5 19

cylinders 30 2π 13.6 11.7 10.9 82clavicles 25 8π 61.9 41.5 157

trachea 23 4π 47.2 33.3 385thrombi 50 2π 39.1 28.0 411

TABLE IIDEAL, RIGID AND NON-RIGID CORRESPONDENCE RESULTS.

C. Rigid versus Non-Rigid Correspondence

The surfaces of each of the six phantom and real populationswere parameterized using the progressive parameterizationtechnique of [14]. The chosen height h of the cylinder foreach of the populations can be found in Table I. From theparameterized surfaces, a 16 × 16 h

2π b-spline representationwas computed for each surface using mp = 10.000 points.This was followed by the construction of the rigid correspon-dence. The number of landmarks to estimate the covariancewas np = 4000 for all populations. Starting from the rigid cor-respondence, the non-rigid correspondence was calculated. Inthis construction, a b-spline regularization factor of αρ = 0.2was used for all populations. Three levels of scale were usedfor the b-spline surface, the reparameterization transformationand for the number of landmarks. On the coarsest scale, a4×4 b-spline surface, together with a 4×4 reparameterizationtransformation, was the covariance matrix was estimated based

12

0.5 1.0 1.5 2.0 2.5number of modes

0.4

0.6

0.8

1.0

1.2

1.4

1.6

cu

mu

lati

ve v

ari

an

ce

x1e+5 Compactness

bsplinerigidideal

1 2 3 4 5number of modes

0.0

0.5

1.0

1.5

2.0

2.5

3.0

3.5

cu

mu

lati

ve v

ari

an

ce

x1e+5 Compactness

bsplinerigidideal

0.5 1.0 1.5 2.0 2.5 3.0 3.5number of modes

1.0

1.5

2.0

2.5

3.0

cu

mu

lati

ve v

ari

an

ce

x1e+5 Compactness

bsplinerigidideal

0.5 1.0 1.5 2.0 2.5number of modes

0.0

0.5

1.0

1.5

2.0

2.5

mean

err

or

(mm

)

Reconstruction

bsplinerigidideal

1 2 3 4 5number of modes

0.0

0.2

0.4

0.6

0.8

1.0

1.2

1.4

mean

err

or

(mm

)

Reconstruction

bsplinerigidideal

0.5 1.0 1.5 2.0 2.5 3.0 3.5number of modes

0.0

0.2

0.4

0.6

0.8

1.0

1.2

1.4

mean

err

or

(mm

)

Reconstruction

bsplinerigidideal

0.5 1.0 1.5 2.0 2.5number of modes

0.0

0.5

1.0

1.5

2.0

2.5

mean

err

or

(mm

)

Generalization

bsplinerigidideal

1 2 3 4 5number of modes

0.0

0.2

0.4

0.6

0.8

1.0

1.2

1.4

1.6

mean

err

or

(mm

)

Generalization

bsplinerigidideal

0.5 1.0 1.5 2.0 2.5 3.0 3.5number of modes

0.0

0.2

0.4

0.6

0.8

1.0

1.2

1.4

1.6

mean

err

or

(mm

)

Generalization

bsplinerigidideal

0.5 1.0 1.5 2.0 2.5number of modes

0.6

0.7

0.8

0.9

1.0

1.1

1.2

mean

err

or

(mm

)

Specificity

bsplinerigidideal

1 2 3 4 5number of modes

0.0

0.5

1.0

1.5

2.0

2.5

3.0

mean

err

or

(mm

)

Specificity

bsplinerigidideal

0.5 1.0 1.5 2.0 2.5 3.0 3.5number of modes

0.30

0.35

0.40

0.45

0.50

0.55

0.60

mean

err

or

(mm

)Specificity

bsplinerigidideal

Fig. 7. The model performance measures, together with the standard error, for the phantom populations. From left to right: beams, disks and cylinderspopulations. From top to bottom: compactness, reconstruction, generalization and specificity measure. In each graph a comparison is made of the idealcorrespondence (green), the rigid correspondence (red) is and the non-rigid correspondence (black). For each model, the first m modes are shown that capture99% of the total variance in the model.

on 250 landmarks. At each new resolution lavel the valuesare increased as detailed in Section VI. All optimizationsproblems were solved with the L-BFGS routine [46]. For thealignment initialization (Eq. 26) and full alignment (Eq. 27)a gradient tolerance of 0.01 was used to determine conver-gence. For the multi-resolution correspondence optimization(Eq. 29) a gradient tolerance of 0.01 × 2−L was used. Thequality of the obtained correspondences was assessed using theabove-mentioned model performance measures. The resultsare shown in Figure 7 for the phantom populations and inFigure 8 for the real populations.

In Figure 7, the ideal (intuitive), the rigid and the non-rigidcorrespondence for the phantom populations are compared.For the disk population, it can be seen that the non-rigid and

the ideal correspondence are of comparable quality. The rigidcorrespondence on the other hand is much worse. This is dueto the fact that the parameterization technique maps the diskpart of the surfaces to a different location in the parameterspace. The non-rigid correspondence improvement, on theother hand, moves the disks to the same part of the parameterspace and this results in an optimal correspondence. See Figure6 for a visualization of the rigid and b-spline optimized model.The improved compactness of the non-rigid over the idealcorrespondence can be attributed to the reduced area that thedisk part of the surface for the non-rigid correspondence takesin the parameter-space. For the beam population, the qualityof the correspondences is comparable, with a slight advantagefor the ideal, since the parameterization technique already

13

5 10 15number of modes

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.0

cu

mu

lati

ve v

ari

an

ce

x1e+5 Compactness

bsplinerigid

2 4 6 8 10number of modes

0.4

0.6

0.8

1.0

1.2

1.4

1.6

1.8

cu

mu

lati

ve v

ari

an

ce

x1e+5 Compactness

bsplinerigid

5 10 15 20 25number of modes

0.5

1.0

1.5

2.0

2.5

3.0

3.5

cu

mu

lati

ve v

ari

an

ce

x1e+5 Compactness

bsplinerigid

5 10 15number of modes

0.2

0.4

0.6

0.8

1.0

1.2

1.4

mean

err

or

(mm

)

Reconstruction

bsplinerigid

2 4 6 8 10number of modes

0.2

0.4

0.6

0.8

1.0

1.2

1.4

1.6

mean

err

or

(mm

)

Reconstruction

bsplinerigid

5 10 15 20 25number of modes

1.0

1.5

2.0

2.5

mean

err

or

(mm

)

Reconstruction

bsplinerigid

5 10 15number of modes

0.6

0.7

0.8

0.9

1.0

1.1

1.2

1.3

1.4

1.5

mean

err

or

(mm

)

Generalization

bsplinerigid

2 4 6 8 10number of modes

0.2

0.4

0.6

0.8

1.0

1.2

1.4

1.6

mean

err

or

(mm

)

Generalization

bsplinerigid

5 10 15 20 25number of modes

1.0

1.5

2.0

2.5

mean

err

or

(mm

)

Generalization

bsplinerigid

5 10 15number of modes

0.95

1.00

1.05

1.10

1.15

1.20

1.25

1.30

mean

err

or

(mm

)

Specificity

bsplinerigid

2 4 6 8 10number of modes

0.7

0.8

0.9

1.0

1.1

mean

err

or

(mm

)

Specificity

bsplinerigid

5 10 15 20 25number of modes

1.7

1.8

1.9

2.0

2.1

2.2

2.3

2.4

2.5

mean

err

or

(mm

)Specificity

bsplinerigid

Fig. 8. Different model performance measures, together with the standard error, for the real data sets. From left to right: clavicles, tracheas and thrombipopulations. From top to bottom: compactness, reconstruction, generalization and specificity measure. In each graph the rigid correspondence (red) is comparedto the non-rigid correspondence (black). For each model, the first m modes are shown that capture 99% of the total variance in the model.

generates a good correspondence. The most notable differenceis the improved specificity of the ideal model, which is dueto the fact that the other models can generate samples withrounded corners. For the cylinders population, it can be seenthat the non-rigid correspondence is an improvement over theideal, which is mainly due to the improved spatial alignmentof the surfaces.

In Figure 6, a visualization of the first few modes ofthe models of the CT-scanned populations can be found. InFigure 8, the model performance measures for the rigid andthe non-rigid correspondence are shown for the CT-scannedpopulations. It is more difficult to analyze these results becausethere is no ideal correspondence available. It can be seenthough that the rigid correspondence generates good modelsfor all three populations and that the non-rigid correspondence

is a significant improvement in most cases. Compared to theclavicle and trachea populations, the approximation errors forthe thrombi population errors are higher. This can be attributedto the large variability that is present in the thrombi population,together with the lower resolution of the ct-scans, that is0.5× 0.5× 0.5mm3 versus 0.5× 0.5× 2.0mm3.

In Table I, the MDL values for the ideal (when available),the rigid, and the non-rigid correspondences can be found forall populations. It can be seen that the b-spline correspondenceoptimization always succeeds in decreasing the MDL-value.Such a decrease indicates that the resulting shape modelbecame less complex. This is most apparent for the diskpopulation: the five modes of the rigid model (MDL-value of17.2) are reduced to a single mode by the b-spline optimization(MDL-value of 4.6). Table I also reports the execution time

14

5 10 15number of modes

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.0

cu

mu

lati

ve v

ari

an

ce

x1e+5 Compactness

0.050.10.20.40.81.6rigid

2 4 6 8 10 12 14 16number of modes

2

3

4

5

6

7

cu

mu

lati

ve v

ari

an

ce

x1e+4 Compactness

3162125250

2 4 6 8 10 12 14 16number of modes

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.0

cu

mu

lati

ve v

ari

an

ce

x1e+5 Compactness

1D-1O-1L3D-1O-1L3D-3O-1L3D-3O-3L

5 10 15number of modes

0.2

0.4

0.6

0.8

1.0

1.2

1.4

mean

err

or

(mm

)

Reconstruction

0.050.10.20.40.81.6rigid

2 4 6 8 10 12 14 16number of modes

0.2

0.4

0.6

0.8

1.0

1.2

1.4

mean

error (

mm

)

Reconstruction

3162125250

2 4 6 8 10 12 14 16number of modes

0.2

0.4

0.6

0.8

1.0

1.2

1.4

mean

error (

mm

)

Reconstruction

1D-1O-1L3D-1O-1L3D-3O-1L3D-3O-3L

5 10 15number of modes

0.6

0.7

0.8

0.9

1.0

1.1

1.2

1.3

1.4

1.5

mean

err

or

(mm

)

Generalization

0.050.10.20.40.81.6rigid

2 4 6 8 10 12 14 16number of modes

0.6

0.7

0.8

0.9

1.0

1.1

1.2

1.3

1.4

1.5

mean

error (

mm

)

Generalization

3162125250

2 4 6 8 10 12 14 16number of modes

0.6

0.7

0.8

0.9

1.0

1.1

1.2

1.3

1.4

1.5

mean

error (

mm

)

Generalization

1D-1O-1L3D-1O-1L3D-3O-1L3D-3O-3L

5 10 15number of modes

0.95

1.00

1.05

1.10

1.15

1.20

1.25

1.30

mean

err

or

(mm

)

Specificity

0.050.10.20.40.81.6rigid

2 4 6 8 10 12 14 16number of modes

0.95

1.00

1.05

1.10

1.15

1.20

1.25

mean

err

or

(mm

)

Specificity

3162125250

2 4 6 8 10 12 14 16number of modes

0.95

1.00

1.05

1.10

1.15

1.20

1.25

mean

err

or

(mm

)Specificity

1D-1O-1L3D-1O-1L3D-3O-1L3D-3O-3L

Fig. 9. Influence of the method parameters on the model performance for the clavicle data set. Left: influence of the b-spline deformation regularizationfactor αρ. Seven models with regularisations from αρ = 0.05 to αρ = 1.6 are compared. The rigid case, i.e. αρ = ∞, is provided as a reference. It canbe observed that models with less regularisation tend to perform better. However, too low values for αρ may result in irregularities. Middle: influence of thenumber of landmarks np used in the estimation of the shape covariance matrix. Four models are shown with number of landmarks from nLp = 31 · 4L−1

to nLp = 250 · 4L−1. It can be seen that the performance is insensitive to this parameter. However, for certain populations, using too few landmarks mayresult in degraded models due to undersampling of highly variable surface regions. Right: influence of the multi-resolution scheme. The scheme encoded asaD-bO-cL uses a levels for the b-spline deformation, b levels for the surface approximation, and c levels for the number of landmarks. It can be observedthat the single resolution scheme 1D-1O-1L generates a degraded model and that the other schemes result in comparable models. The full scheme 3D-3O-3Lis preferred because it is computationally less expensive.

for the construction of the correspondences, once the param-eterizations are obtained. About 10% of the time is taken bythe rigid correspondence construction. The correspondencesfor the phantom populations are constructed more efficientlycompared with the CT-scanned populations. This is becausethe phantom shape models are relatively simple, i.e. theyhave a small number of modes. The construction for the CT-scanned populations takes a couple of hours. Togheter withthe construction of the parameterizations, a correspondencecan easily be established overnight.

D. Influence of Parameters

In this section, the influence of the method parameters onthe resulting correspondence is investigated. Figure 9 showsthe model performance parameters for the clavicle populationwhen using different values for the most important parameters:(a) the b-spline deformation regularization controlled by factorαρ in Eq. (29), (b) the number of landmarks to estimatethe shape covariance matrix, controlled by np in Eq. (14),and (c) the number of scale levels L for the b-spline surfaceapproximation, the b-spline parameterization transform and the

15

number of levelsreparameterization surface landmarks MDL time(h)

1 1 1 42.6 4.53 1 1 41.3 7.53 3 1 41.3 7.53 3 3 41.5 2.5

TABLE IIINFLUENCE OF THE MULTI-LEVEL SCHEME FOR THE CLAVICLE

POPULATION.

landmarks.In the first column of Figure 9, the influence of the b-spline

regularization factor αρ is shown and it can be seen that, asexpected, lower regularization, i.e. smaller αρ, generates bettermodels. However, the regularization can not be lowered toomuch since then irregularities can appear. No irregularitieswere noted with regularization factors αρ ≤ 0.2 for allconsidered populations.

In the second column of Figure 9, it can be seen that theinfluence of the number of landmarks np is negligible forthe clavicle population. It can, however, happen that usingtoo few landmarks results in undersampling of highly variableparts of the surface, which, in turn, will result in a degradedcorrespondence. This was observed for the thrombi population,where nLp = 250 · 4L−1 landmarks generated a significantlybetter correspondence than nLp = 31 · 4L−1 landmarks (resultnot shown). Using less landmarks reduces computation time,but it can also result in a degraded correspondences.

In the last column of Figure 9, the model performancemeasures for different multi-scale schemes are shown andin Table II the corresponding MDL values are listed. Fourdifferent schemes are shown: 1D-1O-1L optimizes the cor-respondence using a single resolution, 3D-1O-1L uses threeresolution levels for the deformations, 3D-3O-1L uses threelevels for the deformations and the surfaces, and 3D-3O-3Lis the full multi-resolution scheme using three levels for thedeformations, the surfaces, and the landmarks. It can be seenthat, for the clavicle, the single scale scheme generates adegraded correspondence and the three other schemes gen-erate comparable correspondences. However, the full 3-scalescheme is considerably faster. For the disk population, theperformance degradation of the single scale method was evenmore notable since, as opposed to the three other schemes, itdid not successfully correspond the disk parts of the surface(result not shown).

VIII. CONCLUSIONS

In this work, the minimum description length approachfor shape modeling was translated to surfaces of cylindricaltopology. The proposed method establishes an alignment anda correspondence for a population of surfaces of cylindricaltopology. It generates a rigid correspondence based on cylin-drical surface parameterizations and an improved correspon-dence using multi-level b-spline reparameterizations. Care wastaken to ensure that the objective functions are differentiablewith respect to the alignment and reparameterization parame-ters and, where necessary, an expression for the gradient was

provided. It was shown that the method produces correspon-dences that agree with the intuitive correspondence and thatthe derived shape models generate small approximation errors.

The cylindrical correspondence method of this paper, to-gether with the spherical correspondence methods from [28]and [29], and the disc-like correspondence method of [30],already cover a wide range of biomedical surfaces. However, itwould be interesting to extend the method to other topologies.The method of this paper can be trivially extended to surfacesof genus-1 topology where the torus can be used as theparametric domain. More complex topologies can be treated bydecomposing the surfaces in a consistent set of discs and tubes.Populations of tubular structures with bifurcations could behandled well with this approach. However, arbitrary complexsurfaces will suffer from the boundary constraints imposed bythe surface decomposition. The development of a method thatcan handle arbitrary topology, without constraints, is a verychallenging problem to be solved in the future.

ACKNOWLEDGEMENTS

This work was financially supported by the Institute forthe Promotion of Innovation through Science and Technologyin Flanders (IWT-Vlaanderen) and the Fund for ScientificResearch (F.W.O)-Flanders, Belgium. The thrombus data setswere kindly provided by Marleen de Bruijne and the Dept. ofVascular Surgery, Univ. Medical Center Utrecht.

REFERENCES

[1] T. F. Cootes, C. J. Taylor, D. H. Cooper, and J. Graham, “Active shapemodels: their training and application,” Computer Vision and ImageUnderstanding, vol. 61, no. 1, pp. 38–59, 1995.

[2] I. L. Dryden and K. Mardia, Statistical Shape Analysis. John WileySons, 1998.

[3] M. Styner, J. A. Lieberman, R. K. McClure, D. R. Weinberger, D. W.Jones, and G. Gerig, “Morphometric analysis of lateral ventricles inschizophrenia and healthy controls regarding genetic and disease-specificfactors,” Proc Natl Acad Sci U S A, vol. 102, no. 13, pp. 4872–4877,Mar 2005.

[4] R. R. Paulsen, R. Larsen, S. Laugesen, C. Nielsen, and B. K. Ersbøll,“Building and testing a statistical shape model of the human ear canal,”in Medical Image Computing and Computer-Assisted Intervention -MICCAI 2002, 5th Int. Conference, Tokyo, Japan,. Springer, 2002.

[5] S. Zachow, H. Lamecker, B. Elsholtz, and M. Stiller, “Reconstruction ofmandibular dysplasia using a statistical 3d shape model,” in Proceedingsof CARS (H. Lemke et al, eds.), Int. Congress Series (1281). Elsevier,2005, pp. 1238–1243.

[6] A. Kelemen, G. Szekely, and G. Gerig, “Elastic model-based segmen-tation of 3-D neuroradiological data sets,” IEEE Trans Med Imaging,vol. 18, no. 10, pp. 828–839, Oct 1999.

[7] T. F. Cootes, A. Hill, C. J. Taylor, and J. Haslam, “Use of active shapemodels for locating structures in medical images,” Image and VisionComputing, vol. 12, no. 6, pp. 355–365, 1994.

[8] M. de Bruijne, B. van Ginneken, M. A. Viergever, and W. J. Niessen,“Interactive segmentation of abdominal aortic aneurysms in cta images,”Medical Image Analysis, vol. 8, no. 2, pp. 127–138, 2004.

[9] M. A. Styner, K. T. Rajamani, L.-P. Nolte, G. Zsemlye, G. Szekely, C. J.Taylor, and R. H. Davies, “Evaluation of 3D correspondence methodsfor model building,” Inf Process Med Imaging, vol. 18, pp. 63–75, Jul2003, evaluation Studies.

[10] F. Bookstein, Morphometric Tools for Landmark Data: Geometry andBiology. Cambridge University Press, 1991.

[11] M. S. Floater, “Parametrization and smooth approximation of surfacetriangulations,” Computer Aided Geometric Design, vol. 14, no. 3, pp.231–250, 1997.

[12] C. Brechbuhler, G. Gerig, and O. Kubler, “Parametrization of closedsurfaces for 3-d shape description,” Comput. Vis. Image Underst.,vol. 61, no. 2, pp. 154–170, 1995.

16

[13] E. Praun and H. Hoppe, “Spherical parametrization and remeshing,”ACM Trans. Graph., vol. 22, no. 3, pp. 340–349, 2003.

[14] T. Huysmans, J. Sijbers, and B. Verdonk, “Parameterization of tubularsurfaces on the cylinder,” The journal of WSCG, vol. 13, no. 3, pp.97–104, 2005.

[15] M. Eck, T. DeRose, T. Duchamp, H. Hoppe, M. Lounsbery, andW. Stuetzle, “Multiresolution analysis of arbitrary meshes,” in SIG-GRAPH ’95: Proceedings of the 22nd annual conference on Computergraphics and interactive techniques. New York, NY, USA: ACM Press,1995, pp. 173–182.

[16] H. Lamecker, T. Lange, and M. Seebaß, “A statistical shape modelfor the liver.” in Medical Image Computing and Computer-AssistedIntervention - MICCAI 2002, 5th International Conference, Tokyo,Japan, September 25-28, 2002, Proceedings, Part II, ser. Lecture Notesin Computer Science, T. Dohi and R. Kikinis, Eds., vol. 2489. Springer,2002, pp. 421–427.

[17] H. Lamecker, M. Seebaß, H.-C. Hege, and P. Deuflhard, “A 3d statisticalshape model of the pelvic bone for segmentation,” in Proceedingsof SPIE - Volume 5370 Medical Imaging 2004: Image Processing,J. Fitzpatrick and M. Sonka, Eds., May 2004, pp. 1341–1351.

[18] R. H. Davies, C. J. Twining, T. F. Cootes, J. C. Waterton, and C. J.Taylor, “A minimum description length approach to statistical shapemodeling,” IEEE Trans Med Imaging, vol. 21, no. 5, pp. 525–537, May2002.

[19] M. K. Hurdal, P. L. Bowers, K. Stephenson, D. W. L. Sumners, K. Rehm,K. Schaper, and D. A. Rottenberg, “Quasi-conformally flat mappingthe human cerebellum,” in MICCAI ’99: Proceedings of the SecondInternational Conference on Medical Image Computing and Computer-Assisted Intervention. London, UK: Springer-Verlag, 1999, pp. 279–286.

[20] C. Gotsman, X. Gu, and A. Sheffer, “Fundamentals of spherical parame-terization for 3d meshes,” in SIGGRAPH ’03: ACM SIGGRAPH 2003Papers. New York, NY, USA: ACM, 2003, pp. 358–363.

[21] X. Gu, Y. Wang, T. Chan, P. Thompson, and S.-T. Yau, “Genuszero surface conformal mapping and its application to brain surfacemapping,” Medical Imaging, IEEE Transactions on, vol. 23, no. 8, pp.949–958, Aug. 2004.

[22] S. Haker, S. Angenent, A. Tannenbaum, and R. Kikinis, “Nondistortingflattening maps and the 3-D visualization of colon CT images,” IEEETrans Med Imaging, vol. 19, no. 7, pp. 665–670, Jul 2000.

[23] W. Hong, X. Gu, F. Qiu, M. Jin, and A. Kaufman, “Conformal virtualcolon flattening,” in Proceedings of the 2006 ACM symposium on Solidand physical modeling. Cardiff, Wales, United Kingdom: ACM, 2006,pp. 85–93.

[24] T. Huysmans and J. Sijbers, “Cylindrical surface parameterization,”submitted.

[25] M. Jin, J. Kim, F. Luo, and X. Gu, “Discrete surface ricci flow,”Visualization and Computer Graphics, IEEE Transactions on, vol. 14,no. 5, pp. 1030–1043, Sept.-Oct. 2008.

[26] A. C. W. Kotcheff and C. J. Taylor, “Automatic construction of eigen-shape models by direct optimization,” Medical Image Analysis, vol. 2,pp. 303–314(12), December 1998.

[27] H. H. Thodberg, “Minimum Description Length Shape and AppearanceModels,” in Information Processing in Medical Imaging, IPMI 2003,2003.

[28] R. H. Davies, C. J. Twining, T. F. Cootes, J. C. Waterton, and C. J. Tay-lor, “3d statistical shape models using direct optimisation of descriptionlength,” in Computer Vision - Eccv 2002 Pt Iii, ser. Lecture Notes inComputer Science, 2002, vol. 2352, pp. 3–20.

[29] T. Heimann, I. Wolf, T. Williams, and H. Meinzer, “3d active shapemodels using gradient descent optimization of description length,” inInformation Processing in Medical Imaging, ser. Lecture Notes inComputer Science, 2005, vol. 3565, pp. 566–577.

[30] P. Horkaew and G.-Z. Yang, “Construction of 3D Dynamic StatisticalDeformable Models for Complex Topological Shapes,” in 7th Interna-tional Conference on Medical Image Computing and Computer-AssistedIntervention (MICCAI 2004), Saint-Malo, France, September 26-29,2004, vol. 3216, September 2004, pp. 217–224.

[31] P. Horkaew and G. Z. Yang, “Optimal deformable surface models for 3dmedical image analysis,” in Information Processing in Medical Imaging,Proceedings, ser. Lecture Notes in Computer Science, 2003, vol. 2732,pp. 13–24.

[32] X. Li, Y. Bao, X. Guo, M. Jin, X. Gu, and H. Qin, “Globally optimal sur-face mapping for surfaces with arbitrary topology,” IEEE Transactionson Visualization and Computer Graphics, vol. 14, no. 4, pp. 805–819,2008.

[33] J. Cates, P. T. Fletcher, M. Styner, M. Shenton, and R. Whitaker, “Shapemodeling and analysis with entropy-based particle systems.” Inf ProcessMed Imaging, vol. 20, no. NIL, pp. 333–45, 2007.

[34] L. Ferrarini, H. Olofsen, W. M. Palm, M. A. van Buchem, J. H. Reiber,and F. Admiraal-Behloul, “Games: Growing and adaptive meshes forfully automatic shape modeling and analysis,” Medical Image Analysis,vol. 11, no. 3, pp. 302 – 314, 2007.

[35] M. de Bruijne, B. van Ginneken, W. Niessen, and M. Viergever, “Adapt-ing Active Shape Models for 3D Segmentation of Tubular Structures inMedical Images,” in Information Processing in Medical Imaging, ser.Lecture Notes in Computer Science, C. Taylor and J. Noble, Eds., vol.2732. Springer, 2003, pp. 136–147.

[36] T. Huysmans, J. Sijbers, F. Vanpoucke, and B. Verdonk, “Improvedshape modeling of tubular objects using cylindrical parameterization,”Lecture Notes in Computer Science, vol. 4091, pp. 84–91, 2006.

[37] R. H. Bartels, J. C. Beatty, and B. A. Barsky, An introduction to splinesfor use in computer graphics & geometric modeling. San Francisco,CA, USA: Morgan Kaufmann Publishers Inc., 1987.

[38] I. Lloyd, N. Trefethen, and D. Bau, Numerical Linear Algebra. 3600University City Science Center, Philadelphia, PA: SIAM, 1997.

[39] I. T. Jolliffe, Principal Component Analysis. Springer, October 2002.[40] A. Ericsson and K. Astrom, “Minimizing the description length using

steepest descent,” in Proc. British Machine Vision Conference, Norwich,United Kingdom, vol. 2, 2003, pp. 93–102.

[41] T. Papadopoulo and M. I. A. Lourakis, “Estimating the jacobian of thesingular value decomposition: Theory and applications.” in ECCV (1),ser. Lecture Notes in Computer Science, D. Vernon, Ed., vol. 1842.Springer, 2000, pp. 554–570.

[42] E. Catmull and J. Clark, “Recursively generated b-spline surfaces onarbitrary topological meshes,” Computer-Aided Design, vol. 10, no. 6,pp. 350–355, November 1978.

[43] P. Horkaew and G.-Z. Yang, “Optimal Deformable Surface Models for3D Medical Image Analysis,” in Information Processing in MedicalImaging, IPMI 2003, 2003, pp. 13–24.

[44] D. Loeckx, “Automated nonrigid intra-patient image registration usingb-splines,” Ph.D. dissertation, K.U.Leuven, Belgium, May 2006.

[45] D. Rueckert, L. I. Sonoda, C. Hayes, D. L. Hill, M. O. Leach, andD. J. Hawkes, “Nonrigid registration using free-form deformations: ap-plication to breast mr images.” IEEE Transactions on Medical Imaging,vol. 18, no. 8, pp. 712–721, August 1999.

[46] D. C. Liu and J. Nocedal, “On the limited memory bfgs method forlarge scale optimization,” Math. Program., vol. 45, no. 3, pp. 503–528,1989.