## Jq

where G* represents convolution with the impulse response filter. As noted in Gramkow and Bro-Nielsen [13], a fast, "demons-based" warping algorithm [106] calculates a flow velocity by regularizing the force field driving the template with a Gaussian filter. Since this filter may be regarded as a separable approximation to the continuum-mechanical filters derived earlier, interest has focused on deriving additional separable (and therefore rapidly applied) filters to capture the deforma-tional behavior of material continua in image registration [55,56]. Variable viscosity fluids create a still more flexible registration model ([71]; X = 0, cf. Eq. (4)):

pV2v(x, t) + ^V(VT • v(x, t)) + (d^/dxj)Vvj(x, t)

In this model, a spatially dependent viscosity field ^(x) controls the pliancy of the deforming image, slowing the deformation of some regions to enable others to register first. By combining a diffusion model with spatially adaptive controls, this method shares some affinities with variable-conductance scale spaces, which have become a powerful approach for multiscale feature detection in medical imaging [105].

### 2.6 Multigrid and Coarse-to-Fine Optimization

Vast numbers of parameters are required to represent complex deformation fields. Optimal values for these parameters must therefore be searched, to find parameter sets that globally minimize the measure of mismatch between the warped image and target. In several robust systems for automated brain image segmentation and labeling, Dengler and Schmidt [32], Bajcsy and Kovacic [6], Collins et al. [20,21], Gee et al. [49,50], and Schormann et al. [97] recover the optimal transformation in a hierarchical multiscale fashion. Both template and target intensity data are smoothed with different sized Gaussian filters, and the registration is performed initially at coarse spatial scales, then finer ones. This accelerates computation and helps to avoid local minima of the mismatch measure. In an alternative approach, the deformation field is expressed as a steadily increasing sum of basis functions, whose coarse-scale parameters are estimated first and whose spatial frequency is increased as the algorithm progresses [1,82]. The widely used Statistical Parametric Mapping [4] and Automated Image Registration algorithms [130] take a similar coarse-to-fine approach. Increasingly complex warping fields are expressed in terms of a 3D cosine basis (SPM) or by using 3D polynomials of increasing order (AIR; polynomials also span the space of continuous deformation fields, by the Stone-Weierstrass theorem). Amit [2] expresses 2D warping fields in terms of a wavelet basis. The small support of high-frequency wavelets allows local adjustments of the warping field in regions where the mismatch is large, without having to alter the field in other areas where the mismatch may already be small. In Miller et al. ([82]; cf. [1]), a stochastic algorithm generates the expansion coefficient set for the deformation field u(x)=Ei,= lto d Er^^2 + j2^1-^W^r M] (8)

in terms of an eigenbasis {ci;j,r} for the linear elasticity operator. Stochastic gradient descent is used to find the optimal warping field parameters according to:

d^,hr (t) = — (1/2)[SH (u(t))/QpiJir]dt + dwhhr (t). (9)

Here H(u(t)) is the combined measure of intensity mismatch and deformation severity, and dw¡,j,r(t) is a Wiener process that allows provisional parameter estimates to jump out of local minima. At the expense of added computation time, stochastic sampling allows globally optimal image matches to be estimated.

### 2.7 Bayesian Registration models

Bayesian statistical models provide a general framework for the presentation of deformable matching methods. They have also been applied to the brain image matching problem [4,49,50,82]. In a Bayesian model, statistical information on the imaging process (the imaging model) is combined with prior information on expected template deformations (the prior model) to allow inferences to be made about the parameters of the deformation field. In fact, all of the intensity-based approaches that combine an intensity mismatch term with a measure of deformation severity can be recast as an inference task in a Bayesian probabilistic framework. The Bayesian maximum a posteriori (MAP) estimator solving the registration problem is the transformation U = argminu(D(u)+ E(u)), which minimizes a combined penalty due to the intensity mismatch D(u) and the elastic potential energy E(u) stored in the deformation field.

### 2.8 Model-Driven Algorithms

The extreme difficulty of matching brain data based on intensity criteria alone led to the development of algorithms driven by anatomical models, which can be extracted from each dataset prior to registration. Anatomic models provide an explicit geometry for individual structures in each scan. Model-driven algorithms can be classified by the type of features that drive them. These features include point, curve, or surface models of anatomic structures. When parameterized consistently (see Section 3), mesh-based models of anatomic systems help guide the mapping of one brain to another. Anatomically driven algorithms guarantee biological as well as computational validity, generating meaningful object-to-object correspondences, especially at the cortex.

### 2.9 Point-Based Matching

The simplest set of anatomic features that can guide the mapping of one brain to another is a set of point landmarks, identified manually [8] or automatically [2,30] in each dataset. A specification of correspondences at point landmarks can be extended to produce a deformation field for the full volume in a variety of ways, each consistent with the displacements assigned at the point landmarks. The ambiguity is resolved by requiring the deformation field to be the one that minimizes a specific regularizing functional [119], which measures the roughness or irregularity of the deformation field, calculated from its spatial derivatives. Different regularizers measure smoothness in different ways, and have different ways of interpolating the deformation from the points into the full 3D volume. All regularizers penalize large values of the derivatives, creating warping functions that take on the values to be interpolated (displacements) at the specified points, but that do not vary erratically elsewhere. For producing 2D warping fields, thin-plate splines [8, 60] are functions that minimize the penalty,

Jthin-piate (u)= / [(811u)2 + 2(812u)2 + (822u)2] dx1 dx2, (10)

JR.2

where u = 82u/8xi8x;-. Other members of the family of multidimensional spline functions [37,80,124] have also been used for brain image warping. Membrane splines [1,49], elastic body splines [30,82], and div-curl splines [102] are warping functions that minimize the following measures of irregularity:

Jmemb(u) = J [(81 u1)2 + (81u2)2 + (82u1)2 + (82u2)2]

^elas (U) = j( Ei=1 to 2 Em to 2 [(^/2)(8^)(8;",)

Just like the continuum-mechanical warps defined earlier, the warping fields generated by splines satisfy partial differential equations of the form Lu(x) = — F(x), where u(x) is fixed at the specified points, F(x) plays the role of a body force term, and L is the biharmonic differential operator V4 for the thin-plate spline, the Laplacian operator V2 for the membrane spline, and the Cauchy-Navier operator pV2 + ( + p)V(VT•) for the elastic body spline.1

Once a type of spline is chosen for warping, a formula can be used that specifies how to interpolate the displacement field from a set of points {x;} to the surrounding 2D plane or 3D volume:

Here pm_ 1 (x) is a polynomial of total degree m — 1, where m is the order of derivative used in the regularizer, and G is a radial basis function or Green's function whose form depends on the type of spline being used [30, 65]. The polynomial coefficients determine components of low-order data variations, such as a global rigid motion between the two brains, which are not penalized by the smoothness functional [127]. The coefficients ci in this formula are found by evaluating this equation at the points xi and solving the resulting linear system of equations. Various choices of radial basis functions G(r), for r =||(x — xi)||, are commonly used for medical image matching and for multivariate data interpolation generally [8, 94, 111]. Their behavior is also well understood from their wide use in neural networks and pattern classifiers [92]. Choices of r2 ln r and r correspond to the thin-plate spline in two and three dimensions, with r3 for the 3D volume spline [30], the 3x3 matrix [ar21 — 3xxr]r for the 3D elastic body spline [30], and (r2 + c2)a(0<a< 1) and ln(r2 + c2)1/2(c2 > 1) for multi-quadric and shifted-log interpolants [42,63,94]. Gaussian radial functions G(r) = exp(—r2/2a2) belong to the same family, generating warping fields that minimize the following nonintuitive irregularity measure [89]:

1 Note that spline-based warping functions can be defined either (1) as the variational minimizer of an irregularity measure or (2) as solutions to related PDEs. But, as pointed out by Arad [3], strictly one looks for solutions to each problem in different function spaces (the Sobolev space H2 (fl) for minimizing the regularizing integral — which is also known as a Sobolev semi-norm — and continuous function spaces such as C(flc) n C4(fl) for solving PDEs such as V4u(x) = 0, x e fl).

These Gaussian functions perform poorly unless a is chosen carefully ([127]; Franke [42] suggests using a = 1.008d/n, where d is the diameter of the point set and n is the number of points). The best choice of radial basis functions depends on the objective. For modeling deformation of real biological tissue, elastic body splines outperform thin-plate splines [30]. Thin-plate splines, however, are advantageous in computing shape statistics [8]. Use of radial basis functions in medical uk(x) = ^^ image registration algorithms has also been extended to neural network implementations [29].

functions, or convolution filters, in the continuum-mechanical matching approach [13, 65]. They are also directly analogous to Watson-Nadaraya kernel estimators, or Parzen windows, in nonparametric regression methods ([87]; Fig. 2). As before, the k deformation field components are the output values of the neural net:

2.10 Neural Network Approaches

Neural networks are widely used in pattern recognition and computational learning [92] and provide ingenious solution for both landmark-driven and automated image matching. Two of the most common network types are multilayered perceptrons (MLPs) and radial basis function (RBF) networks. RBF neural nets, in particular, use correspondences at known landmarks as a training set to learn a multivariate function mapping positions in the image (input) to the desired displacement field at that point (output). Intriguingly, the hidden units in the neural net are directly analogous to Green's

Here the Gj are N separate hidden unit neurons with receptive fields centered at x{, J2

am^^mm is a polynomial whose terms are hidden units and whose coefficients am are also learned from the training set, and wik are synaptic weights (Fig. 2). The synaptic weights are determined by solving a linear system obtained by substituting the training data into this equation. If landmarks are available to constrain the mapping, the function centers xi may be initialized at the landmark positions, although fewer hidden units are desirable if the network output is to be accurate far from the landmarks. In an innovation [29], a 3D brain image matching neural net was developed that eliminates the need for landmark selection. Network weights (the coordinate transformation parameters) and the RBF center locations are successively tuned to optimize an intensity-based functional (normalized correlation) that

FIGURE 2 Automated image matching with neural networks. Neural networks can be used to compute the deformation field required to warp one image into the shape of another. In one approach ([30], (a)), each of the three deformation vector components, uk(x), is the output of the neural net when the position in the image to be deformed, x, is input to the net. In calculating the output, the numerical outputs of the hidden units (Gj, nm) are weighted using synaptic weights, wa, different for each output component. If landmarks in both images are available to guide the mapping, the weights are found by solving a linear system. Otherwise, the weights can be tuned so that a measure of similarity between the deforming image and the target image is optimized. (b) When the number of landmarks driving the deformation mapping is very large [115], a linear system can be set up, based on a model of linear elasticity, which is solved by successive overrelaxation methods. However, there is a strong mathematical connection between the Green's function (impulse response) of the deformation operator (b), hidden units in the neural net (a), and kernels used in nonparametric regression approaches (see main text). (Diagram (a) is adapted from [30]).

FIGURE 2 Automated image matching with neural networks. Neural networks can be used to compute the deformation field required to warp one image into the shape of another. In one approach ([30], (a)), each of the three deformation vector components, uk(x), is the output of the neural net when the position in the image to be deformed, x, is input to the net. In calculating the output, the numerical outputs of the hidden units (Gj, nm) are weighted using synaptic weights, wa, different for each output component. If landmarks in both images are available to guide the mapping, the weights are found by solving a linear system. Otherwise, the weights can be tuned so that a measure of similarity between the deforming image and the target image is optimized. (b) When the number of landmarks driving the deformation mapping is very large [115], a linear system can be set up, based on a model of linear elasticity, which is solved by successive overrelaxation methods. However, there is a strong mathematical connection between the Green's function (impulse response) of the deformation operator (b), hidden units in the neural net (a), and kernels used in nonparametric regression approaches (see main text). (Diagram (a) is adapted from [30]).

measures the quality of the match. Hidden units were initially randomly placed across the image and the network was trained (i.e., the parameters of the warping field were determined) by evaluating the gradient of the normalized correlation with respect to the network parameters, and optimizing their values by gradient descent. Results matching 3D brain image pairs were impressive [29]. For further discussion of the close relationship among continuum-mechanical PDEs, statistical regression, and neural nets, see Ripley et al. [92].

### 2.11 Curve-Based Approaches

When constructing a warping field for matching two brain images, greater anatomical accuracy can be achieved by using curves as anatomic driving features. In fact, many investigators of point-based landmark matching have included orientation attributes or directional information at landmarks to further constrain the transformation [9,77]. In two dimensions, matching of entire planar curved boundaries is ideal for correcting deformations in histologic tissue ([79,97]; see Fig. 3). Approaches using sulcal lines to drive a 3D volumetric warp are under active investigation in Macaque [65] and human MR data [7, 31, 48, 72]. Curve-based sulcal models can be combined with intensity-based measures to assist in matching cortical regions [22].

Declerck et al. [31] used crest lines to drive a volume transformation. Crest lines [83] are curved loci on a surface that satisfy the following geometric criterion: The largest principal curvature must be locally maximal in the associated principal direction. After an MR dataset is thresholded to segment the cortex and ventricles, crest lines on these surfaces are defined. They are matched with their counterparts in a target brain using (1) an iterative closest point algorithm, which finds candidate lines for matching, and (2) topological criteria to enforce one-to-one matching of curves and to ensure that their internal points are matched in a consistent, serial order. The 3D deformation field, expressed as a 3D tensor product of B-spline basis functions, is obtained by minimizing a regularizing term (as in point-based approaches) and a landmark mismatch term. This mismatch term limits the impact of spurious matches by tolerating some separation between curves that the algorithm decides to match.

### 2.12 Automated Matching

Identifying the subset of curved features that have a consistent topology across subjects (and are therefore appropriate to match) is extremely challenging. The problem is, however, easier than in the general intensity-matching case, as the parametric forms of the curves have a differentiable structure. The inherent structure in the model allows additional geometric features (torsion, curvature, and local Frenet frames) to be included in the matching function, to favor correct pairing [54,70]. Sulcal curvature patterns are only weakly conserved across subjects, and even this consistency is restricted to specific anatomic regions [86,113]. Nevertheless, to help guide the automated matching of curves and surfaces in anatomic data, statistical priors have been defined, for different types of sulci, to describe their expected curvature [67,75], torsion [61,67], and stereotaxic position [109,112,115]. In an alternative approach based on Markov random fields, Mangin

FIGURE 3 Curve-driven warping in a histology application. Curve-driven warping algorithms can recover and compensate for patterns of tissue change which occur in post mortem histologic experiments. A brain section ((a), left panel), gridded to produce tissue elements for biochemical assays, is reconfigured ((a), right) into its original position in the cryosection blockface ([79], algorithm from [111], modified to produce a 2D elastic warp). The complexity of the required deformation vector field in a small tissue region ( magnified vector map, (b)) demonstrates that very flexible, high-dimensional transformations are essential [87, 111]. These data can also be projected, using additional warping algorithms, onto in vivo MRI and coregistered PET data from the same subject for digital correlation and analysis [79]. See also Plate 82.

FIGURE 3 Curve-driven warping in a histology application. Curve-driven warping algorithms can recover and compensate for patterns of tissue change which occur in post mortem histologic experiments. A brain section ((a), left panel), gridded to produce tissue elements for biochemical assays, is reconfigured ((a), right) into its original position in the cryosection blockface ([79], algorithm from [111], modified to produce a 2D elastic warp). The complexity of the required deformation vector field in a small tissue region ( magnified vector map, (b)) demonstrates that very flexible, high-dimensional transformations are essential [87, 111]. These data can also be projected, using additional warping algorithms, onto in vivo MRI and coregistered PET data from the same subject for digital correlation and analysis [79]. See also Plate 82.

## Post a comment