2.1 Model-Driven and Intensity-Driven

A wide variety of 3D image warping algorithms have been designed to handle neuroanatomic data (Table 1). Model-driven algorithms first build explicit geometric models, < representing separate, identifiable anatomic elements in ; each of the scans to be matched. These anatomical systems typically include functionally important surfaces

TABLE 1 Warping algorithms for 3D nonlinear deformation of brain data Algorithm type Author General features

[26,34,84,103,111], curves [48,83,101], and point landmarks [1,8]. Anatomical elements are parameterized and matched with their counterparts in the target scan, and their correspondences guide the volumetric transformation of one brain to another. In our own warping algorithms (Section 3; [111,115]), higher-level structural information guides the mapping of one brain onto another, and a hierarchy of curve-to-curve and surface-to-surface mappings is set up, guaranteeing the biological validity of the resulting transform. The algorithms exploit anatomical information to match cortical regions, so that networks of sulci and gyri are

I. Intensity-Based

A. Elastic Broit Connection developed between 3D image and a mechanical system

Bajcsy/Kovacic Multiscale/multiresolution approach; matching based on normalized cross-correlation in a 3D Hermite basis

Dengler/Schmidt "Dynamic pyramid'Voptical flow approach; matching performed on the sign of pixels in a multiresolution Laplacian pyramid; used for template-driven segmentation Miller/Christensen/ MASPAR (massively parallel) implementation; Bayesian (probabilistic) pattern-theoretic framework; warp

Grenander represented using eigenfunctions of elastic operator and solved in a stochastic PDE framework

Gee Bayesian (probabilistic) framework; uses finite elements; matches curvature and edge features as well as intensities and tissue classes Schormann/Zilles Fast elastic matching for histologic applications

B. Hyperelastic Rabbitt Allows large deformation and temporal evolution of the anatomic template

C. Spline-Based Meyer/Kim Mutual information measures pattern similarity; allows cross-modality warping

D. Fluid Christensen/Miller/ MASPAR (massively parallel) implementation; allows large deformations; deformation velocity tracked in

Grenander an Eulerian reference system; topology of the volumetric mapping investigated (related work by

Freeborough, Lester)

Bro-Nielsen/Gramkow Fast filters used to implement continuum-mechanical deformations

E. Other Collins Maximizes normalized cross-correlation of intensity and edge features; used for automated image segmentation and labeling, and probabilistic atlas construction Thirion Fast algorithm; "demon"-based diffusion method motivated by thermodynamics

Woods Automated polynomial-based image registration

Ashburner/Friston Bayesian (probabilistic) framework; represents warp using 3D trigonometric basis; part of Statistical

Amit Wavelet, and other, bases used to express deformation fields in a variational framework, graphical templates developed for automated feature detection

II. Model-Based Constraints

A. Points Bookstein Thin-plate splines; used to investigate biological shape variation

Davis Landmarks used to drive 3D volume splines, elastic splines; neural network implementation

Müller/Ruprecht Explores 3D interpolation schemes and radial basis functions for warping

B. Curves Subsol/Declerck Pre-extracted 'crest lines' constrain the warp; used for automated atlas construction

Ge; Banerjee; Collins 3D sulcal lines used for cortical surface matching

C. Surfaces Downs Normalizes cortical convex hulls

Terzopoulos Snakes used for deformable curve and surface segmentation

Szeliski/Lavallee Warp rapidly computed on an adaptive octree-spline grid

Davatzikos Deformable surface models drive a 3D elastic registration; curvature maps used to constrain surface-to-

surface matching (related work by Gabrani/Tretiak) Thompson/Toga Connected systems of deep surface meshes and deformable surfaces used to drive the warp; tensor-based matching of cortical patterns for pathology detection; used to measure brain growth and to create disease-specific brain atlases

Drury/Van Essen; Dale/ Cortical flattening and parameterization algorithms (related work by Schwartz/Merker, MacDonald/ Sereno/Fischl Evans)

Miller/Drury/Van Essen Matching of cortical flat maps

D. Other Joshi Generalized Dirichlet problem formulated for point, curve, surface and volumetric matching individually matched. These strategies are discussed in Section 3.

Model-driven approaches contrast with intensity-driven approaches. Intensity-driven approaches aim to match regional intensity patterns in each scan based on mathematical or statistical criteria. Typically, they define a mathematical measure of intensity similarity between the deforming scan and the target. Measures of intensity similarity can include squared differences in pixel intensities [4,15,128,130], regional correlation [6,21], or mutual information [69]. Mutual information has proved to be an excellent similarity measure for cross-modality registrations, since it assumes only that the statistical dependence of the voxel intensities is maximal when the images are geometrically aligned. The intensity similarity measure, combined with a measure of the structural integrity of the deforming scan, is optimized by adjusting parameters of the deformation field. Algorithms based on intensity patterns alone essentially bypass information on the internal topology of the brain. Matching of neuroanatomic data in the absence of higher-level structural information presents an extremely difficult pattern recognition problem. Future hybrid approaches, based on a combination of modelbased and densitometric criteria, are likely to benefit from the advantages of each strategy.

A key insight, which spurred the development of intensity-based warping algorithms, was the connection of the image data with a physically deforming system in three dimensions [12]. Physical continuum models (Fig. 1) consider the deforming image to be embedded in a three-dimensional deformable medium, which can be either an elastic material or a viscous fluid. The medium is subjected to certain distributed internal forces, which reconfigure the medium and eventually lead the image to match the target. These forces can be based mathematically on the local intensity patterns in the datasets, with local forces designed to match image regions of similar intensity.

In elastic media, the displacement field u(x) resulting from internal deformation forces F(x) (called body forces) obeys the Navier-Stokes equilibrium equations for linear elasticity:

Here R is a discrete lattice representation of the scan to be transformed, VT • u(x) = ^ 0uj/0xj is the cubical dilation of the medium, V2 is the Laplacian operator, and Lame's coefficients X and ^ refer to the elastic properties of the medium. Body forces, designed to match regions in each dataset with high intensity similarity, can be derived from the gradient of a local correlation function. In Bajcsy and Kovacic

[6], intensity neighborhoods to be correlated in each scan were projected onto a truncated 3D Hermite polynomial basis to enhance the response of edge features and accelerate computation. More complex local operators can also be designed to identify candidates for regional matches in the target image [2]. With proper boundary conditions the elasticity equilibrium equations can be solved numerically by finite difference, finite element, or spectral methods. If U: x^-x + u(x), then U(R) is the final image, warped into register with the target scan. This elasticity-based warping scheme was introduced by Broit [12]. It was subsequently extended to a multiresolution/multigrid algorithm by Bajcsy and Kovacic [6], and to a finite element implementation by Gee et al. [49].

More recently, Christensen et al. [15-17] proposed a viscous-fluid based warping transform, motivated by capturing nonlinear topological behavior and large image deformations. Designed to operate sequentially, this transform is actually a series of three algorithms that adjust successively finer features of the local anatomy until the transformed image matches the target scan. The optimal deformation field is defined as the one that maximizes a global intensity similarity function (defined on the deformed template and the target), while satisfying additional continuum-mechanical constraints that guarantee the topological integrity of the deformed template.

The transformation proposed by Christensen et al. [17] is conducted in three successive stages. Stage 1 requires a sparse manual specification of the displacement field by isolating several corresponding landmarks in both 3D scans. The minimum energy configuration of the template compatible with this initial assignment is formalized as a Dirichlet problem [65]. For a system ofpoint landmarks, the associated Fredholm integral equation reduces to a linear system whose solution expresses the stage 1 deformation field in terms of the self-adjoint linear operator describing the mechanics of the deforming system. A second step expresses the residual deformation in terms of an approximation series of eigenfunc-tions of the linear elastic operator. Basis coefficients are determined by gradient descent on a cost functional (2) that penalizes squared intensity mismatch between the deforming template T(x — u(x, t)) and target S(x):

C(T(x), S(x), u) = (1/2) / |T(x — u(x, t)) — S(x)\2dx. (2) Jo

Finally, a third, viscous deformation stage allows large-distance nonlinear fluid evolution of the neuroanatomic template. With the introduction of concepts such as deformation velocity and an Eulerian reference frame, the energetics of the deformed medium are hypothesized to be relaxed in a highly viscous fluid. The driving force, which deforms the anatomic template, is critical for successful registration. It is

FIGURE 1 Continuum-mechanical warping, (a) The complex transformation required to reconfigure one brain into the shape of another can be determined using continuum-mechanical models, which describe how real physical materials deform, In this illustration, two line elements embedded in a linearly elastic 3D block (lower left) are slightly perturbed (arrows), and the goal is to find how the rest of the material deforms in response to this displacement. The Cauchy-Navier equations (shown in discrete form, top) are solved to determine the values of the displacement field vectors, u(x), throughout the 3D volume. (b) Lame Elasticity Coefficients. Different choices of elasticity coefficients, X and in the Cauchy-Navier equations (shown in continuous form, top) result in different deformations, even if the applied internal displacements are the same. In histologic applications where an elastic tissue deformation is estimated, values of the elasticity coefficients can be chosen which limit the amount of curl (lower right) in the deformation field. Stiffer material models ( top left) may better reflect the deformational behavior of tissue during histologic staining procedures. Note: For visualization purposes, and to emphasize differences in deformation patterns, the magnitudes of the displacement vector fields shown in this figure have been multiplied by a factor of 10. The Cauchy-Navier equations, derived using an assumption of small displacements, are valid only when the magnitude of the deformation field is small.

defined as the variation of the cost functional with respect to the displacement field:

F((x, u(x, t)) = - (T(x - u(x, t)) - S(x))VT\x_u{x t) (3) ^V2v(x, t) + ( + ^)V(VT • v(x, t)) + F(x, u(x, t)) = 0 (4) 8u(x, t)/8t = v(x, t) - Vu(x, t)v(x, t). (5)

The deformation velocity (4) is governed by the creeping flow momentum equation for a Newtonian fluid and the conventional displacement field in a Lagrangian reference system (5) is connected to an Eulerian velocity field by the relation of material differentiation. Experimental results were excellent [17].

Linkage of continuum-mechanical models with a 3D intensity matching optimization problem results in an extremely computationally intensive algorithm. Registration of two 1283 MR volumes took 9.5 and 13 hours for elastic and fluid transforms, respectively, on a 128 x 64 DECmpp1200Sx/ Model 200 MASPAR (massively parallel mesh-connected supercomputer). This spurred work to modify the algorithm for use on standard single-processor workstations [13,43,106].

Both elastic and fluid algorithms contain core systems of up to 0.1 billion simultaneous partial differential equations (1) and (4), requiring many iterations of successive overrelaxation to find their solution. To avoid this need to pass a filter many times over the vector arrays, Bro-Nielsen and Gramkow [13] developed a convolution filter to solve the system of partial differential equations in a single pass over the data. This speeds up the core step in the registration procedure by a factor of 1000. Since the behavior of the mechanical system is governed by the Navier-Stokes differential operator L = ^V2 + (X + ^)V(VT•), the eigenfunctions of this operator [17] were used to derive a Green's function solution u*(x) = G(x) to the impulse response equation Lu*(x) = 5(x — x0). The solution to the full PDE Lu (x) = — F (x) was approximated as a rapid filtering operation on the 3D arrays representing the components of the body force, u(x) = — / G(x — r) • F(r)dr =— (G * F)(x), (6)

Was this article helpful?

## Post a comment