Footnote 3. Deformable surface models often have difficulty in extracting surfaces of structures with extremely complex geometry, such as the ventricles and caudate, due to their irregular or 3D spiral nature, and modifying their behavior is the topic of active research (Cohen and Cohen, 1992; Davatzikos and Prince, 1996; Xu and Prince, 1996). Interestingly, many 3D image warping algorithms are capable of diffeomorphically transforming an irregular or spiral structure into a more regular shape. This more regular shape can then be used as the target for automated surface extraction. The resulting regular surface model can then be back-transformed (using the inverse of the transformation used to warp the underlying image) to reflect its actual irregular geometry. Alternatively, by applying the inverse of the derived 3D deformation field to a regular 3D lattice ruled over the warped 3D image, a reparameterization of the original image with a 3D curvilinear computational mesh is obtained. Since the Riemannian metric tensor of this reparameterization is known, the PDE governing surface extraction can be rewritten in covariant notation (see later, Section 3) and run directly on the original image, using accessory 3D arrays of Christoffel symbols and metric tensor coefficients, obtained from the 3D warping algorithm, to correct the forms of the differential operators. We are currently testing this locally-adaptive modification of the core surface extraction algorithm. In particular, the solution of the surface extraction equation on a 3D curvilinear mesh biases the dynamics of the deforming surface in favor of the expected geometry of the surface to be extracted, without the need to resample the image in which the surface evolves.

Colorized RGB maps of the Cortical Parameter Space. Because the cortical model is arrived at by the deformation of a spherical mesh, any point on the cortical surface (Fig. 8(a)) must map to exactly one point on the sphere (Fig. 8(b)), and vice versa. Each cortical surface is parameterized with an invertible mapping D_p,D_q: (r,s)-> (x,y,z), so sulcal curves and landmarks in the folded brain surface can be reidentified in the spherical map (cf. Sereno et al., 1996 for a similar approach; see also Footnote 4 and Fig. 8(b)). To retain relevant 3D information, cortical surface point position vectors in 3D stereotaxic space are color-coded at 16 bits per channel, to form an image of the parameter space in RGB color image format [Fig. 8(b)/(c)]. To find good matches between cortical regions in different subjects [Fig. 8(a)/(d)], we first derive a spherical map for each respective surface model [Fig. 8(b)/(c)] and then perform the matching process in the spherical parametric space, W . The parameter shift function u(r):W -> W , is given by the solution F_pq:r->r+u(r) to a curve-driven warp in the biperiodic parametric space W =[0,2.pi)x[0,pi) of the external cortex (Fig. 9; cf. Davatzikos, 1996; Drury et al., 1996a; Van Essen et al., 1997).

Fig. 9. Matching Gyral Patterns from One Subject to Another, to Measure Cortical Variability.

Footnote 4: Color-coded spherical maps store geometric and topological information about the sulcal pattern in a 2D array structure which simplifies feature localization and comparison. Cortical models can also be flattened into a fully 2D planar representation (Van Essen and Maunsell, 1990; Schwartz and Merker, 1986; Carman et al., 1995; Drury et al., 1996; Van Essen et al., 1997). Here, the challenging issue of creating planar 2D maps which optimally preserve point-to-point distances is avoided, because the 16-bit color-encoding of the spherical map enables measures of in-surface distance between points on the cortex to computed from the 3D cortical model directly, rather than from a derived or flattened map of the cortical parameter space.

Cortical Curvature. At first sight it seems that the previous approach, using ordinary differential operators to constrain the mapping, can be applied directly. For points r=(r,s) in the parameter space, a system of simultaneous partial differential equations can be written for u(r) using:

L(u(r)) + F(r-u(r)) = 0, for all r in W ,

u(r) = u₀(r) , for all r in M₀ or M₁.

Here M₀, M₁ are sets of points and (sulcal and gyral) curves where vectors u(r) = u₀(r) matching regions of one anatomy with their counterparts in the other are known, and L and F are 2-dimensional equivalents of the differential operators and body forces defined above. Unfortunately, the recovered solution x->D_q(F_pq(D_p^-1(x))) will in general be prone to variations in the metric tensors g_jk(r_p) and g_jk(r_q) of the mappings D_p and D_q. Since the cortex is not a developable surface (Davatzikos, 1996), it cannot be given a parameterization whose metric tensor is uniform (see Footnote 5). As in fluid dynamics or general relativity applications, the intrinsic curvature of the solution domain should be taken into account when computing flow vector fields in the cortical parameter space, and mapping one mesh surface onto another.

Footnote 5. A tensor-based model of the cortex was introduced for biometric applications in a theoretical paper by Bookstein (Bookstein, 1995). In (Thompson and Toga, 1998), we introduce a covariant formalism as a means to directly incorporate metric tensor information into the regularization equations governing the cortical matching process.

Covariant Formalism. To counteract this problem, we developed a covariant formalism (Thompson and Toga, 1998) which makes cortical mappings independent of how each cortical model is parameterized. This mathematical strategy was introduced by Einstein (1914) to allow the solution of physical field equations defined by elliptic operators on manifolds with intrinsic curvature. Similarly, the problem of deforming one cortex onto another involves solving a similar system of elliptic partial differential equations (Drury et al., 1996; Davatzikos, 1996; Thompson and Toga, 1998), defined on an intrinsically curved computational mesh (in the shape of the cortex). In the covariant formalism, the forms of the differential operators governing the mapping of one cortex to another are modified in an adaptive way which reflects the changes in the underlying metric tensor of the surface parameterizations.

Summary of Covariant Approach to Cortical Matching. Cortical surfaces are matched as follows. We first establish the cortical parameterization, in each subject, as the solution of a time-dependent partial differential equation (PDE) with a spherical computational mesh (Thompson and Toga, 1996; Davatzikos, 1996). This procedure sets up an invertible parameterization of each surface in deformable spherical coordinates (Fig. 9), and allows direct computation of the metric tensors g_jk(r_p) and g_jk(r_q) of the mappings. The solution to this PDE defines a Riemannian manifold (Bookstein, 1995). In contrast to prior approaches, this Riemannian manifold is then not flattened (as in Drury et al., 1996; Van Essen et al., 1997), but is used directly as a computational mesh on which a second PDE is defined (see Fig. 10). The second PDE matches sulcal networks from subject to subject. Dependencies between the metric tensors of the underlying surface parameterizations and the matching field itself are eliminated through the use of generalized coordinates and Christoffel symbols (Thompson and Toga, 1998). In the PDE formulation, we replace L by the covariant differential operator L^#. In L^#, all L's partial derivatives are replaced with covariant derivatives. These covariant derivatives are defined with respect to the metric tensor of the surface domain where calculations are performed.

Fig. 10. Covariant Tensor Approach to Cortical Matching.

The covariant derivative of a (contravariant) vector field, uⁱ(x), is defined as (Schutz, 1990):

uⁱ_,k = du^j/dx^k + G ^j_ik uⁱ (20)

where d represents partial derivative, and the G ^j_ik are Christoffel symbols of the second kind. This expression involves not only the rate of change of the vector field itself, as we move along the cortical model, but also the rate of change of the local basis, which itself varies due to the intrinsic curvature of the cortex (cf. Joshi et al., 1995b). On a surface with no intrinsic curvature, the extra terms (Christoffel symbols) vanish. The Christoffel symbols are expressed in terms of derivatives of components of the metric tensor g_jk(x), which can be calculated directly from the cortical model:

G ⁱ_jk = (1/2) g^il (dg_lj/dx^k+dg_lk/dx^j-dg_jk/dxⁱ ). (21)

Scalar, vector and tensor quantities, in addition to the Christoffel symbols required to implement the diffusion operators on a curved manifold are evaluated by finite differences. These correction terms are then used in the solution of the Dirichlet problem for matching one cortex with another (see Footnote 6).

Footnote 6. Several factors complicate the direct application of covariant regularization equations to cortical anatomic data. In fluid mechanics and general relativity applications, usually a single metric tensor field relates the computation mesh to the physical domain of the data, and the solution to the problem is then back-transformed from the computational mesh to the physical domain. In this application however, different metric tensors g_jk(r_p) and g_jk(r_q) relate (1) the physical domain of the input data to the computation mesh (via mapping D_p^-1), and (2) the solution on the computation mesh to the output data domain (via mapping D_q). To address this problem (Fig. 10), the PDE L^##u(r_q)=-F is solved first, to find a flow field T_q:r->r-u(r) on the target spherical map with anatomically-driven boundary conditions u(r_q) = u₀(r_q), for all r_q in M₀ or M₁. Here L^## is the covariant adjustment of the differential operator L with respect to the tensor field g_jk(r_q) induced by D_q. Next, the PDE L^#u(r_p)=-F is solved, to find a reparameterization T_p:r->r-u(r) of the initial spherical map with boundary conditions u(r_p) = 0, for all r_p in M₀ or M₁. Here L^# is the covariant adjustment of L with respect to the tensor field g_jk(r_p) induced by D_p. The full cortical matching field (Fig. 9) is then defined as x->D_q(F_pq(D_p^-1(x))) with F_pq = (T_q) ^-1o (T_p) ^-1_.

Advantages. Using this method, surface matching can be driven by anatomically significant surface features, and the mappings are independent of the chosen parameterization for the surfaces being matched. High spatial accuracy of the match is guaranteed in regions of functional significance or structural complexity, such as sulcal curves and cortical landmarks (Fig. 8). Consequently, the transformation of one cortical surface model onto another is parameterized by one translation vector for each mesh point in the surface model, or 3x65536 ~ 0.2 million parameters. This high-dimensional parameterization of the transformation is required to accommodate fine anatomical variations (cf. Christensen et al., 1995).

This concludes our discussion of warping algorithms. Next, we see how warping algorithms can be used as image analysis tools to support pathology detection in individual subjects or groups.

Related Publications