Mapping techniques for aligning sulci across multiple brains^☆

3.1. Surface inflation

Surface inflation has been previously utilized for visualization purposes (Fischl et al., 1999; Sereno et al., 1996). In our work, surface inflation is used as an initial step, prior to conformal mapping. By using surface inflation prior to conformal mapping, we markedly reduce the dependency on local mesh structure during conformal mapping. To obtain a more regular mesh structure, the surface mesh is smoothed using a relaxation operator (Drury et al., 1996; Timsari and Leahy, 2000); specifically, the vertices are iteratively repositioned according to

where V_i is the position of the ith vertex, t is the iteration number, λ∈[0,1] is a smoothing parameter, and V̄_i is the average vertex position, defined as follows. Let N_i be the set of all triangles that contains vertex V_i, and let C_j and B_j be the center and area, respectively, of triangle j. Then

Surface inflation is a parametric surface deformation method. Our goal is to stop relaxation at the point of a partially flattened surface that should be comparable across individuals. Accordingly, we must define a stopping criterion based on the geometry of the desired final surface. The relaxation operator defined in Eqs. (1) and (2) smooths the initial surface, and mean curvature is one of the geometric features to quantify the amount of smoothing. During relaxation, however, the surface decreases in size. A global measure that quantifies the global shape of the surface and allows comparison between cortical surfaces from different individuals should be scale invariant to capture information about the shape of the surface. Our stopping criterion is based on the L₂ norm of mean curvature, ||H||₂ (Smith et al., 2000), which is defined as

where H is the mean curvature. This is a global measure of extrinsic geometry, and is minimum for a sphere.

The measure ||H||₂ is approximated by discrete differential geometry operators (Meyer et al., 2002). Let κ̄_i = H_in_i be the mean curvature normal at vertex V_i, where H_i is the mean curvature, and n_i is the surface normal. The mean curvature normal κ̄_i can be expressed as

where α_ij and γ_ij are the two angles opposite to the edge in the two triangles sharing the edge (V_i,V_j), as shown in Fig. 4(a), and A_Vi is the surface area of the 1-ring neighborhood (N_i) of vertex V_i with piecewise linear boundaries crossing the triangle edges at their midpoints with right angle, as illustrated in Fig. 4(b).

We stop the relaxation when ||H||₂ ≤ β for a pre-selected β. The left column of Fig. 5 shows an initial cortical surface and several stages during relaxation corresponding to stopping criteria of β=3.25, β=2.5 and β=2.0. The corresponding distributions of triangle angles in the meshes are shown to the right of the figures. A regular mesh should have angles more nearly equal to 60°, which is apparently approximately satisfied by all three partially flattened maps. In order to preserve the most prominent anatomical details representing the major sulci, we have selected β=3.25 as a stopping criterion in this paper [which corresponds to Fig. 5(b)]. In future work, we expect to optimize this selection rather than using this arbitrary choice.

3.2. Conformal mapping

Angenent et al. (1999) proposed a way to construct a conformal flattening equivalence as the solution of a second-order elliptic partial differential equation on the 3D surface to be flattened. Solution of this partial differential equation represents an analytic procedure to map a cortical surface to the complex plane; for triangulated surfaces, a numerical procedure involving a finite element approximation can be used. We have implemented and augmented this basic approach in order to make the resulting maps more standard and have minimum metric distortion (Tosun and Prince, 2001), as we now describe.

The left and right cortical hemispheres are identified by defining a cut around the corpus callosum, as shown in Fig. 6(a) (the red path). This is done automatically using knowledge of the locations of the anterior and posterior commissures (the AC and PC points), which are previously identified for Talairach registration (Talairach and Tournoux, 1988). The line through the AC and PC points is aligned with a vertical line parallel to the y-axis (points from the tip of the patient’s nose to the back of his/her head) during resampling of the MR image volume to obtain isotropic voxels. This line defines the location of the “mid-sagittal” vertical plane. The vertices of triangles cut by the mid-sagittal plane form the hemispherical cut path around the corpus callosum. A fixed north pole is required for conformal mapping using the method of (Angenent et al., 1999). The north pole point is chosen automatically on the top of the corpus callosum, in a triangle adjacent to the hemispherical cut path and equidistant to the AC and PC points, as shown in Fig. 6(a) (yellow point). The entire surface (both hemispheres) are then mapped to the complex plane using the technique described in (Angenent et al., 1999). Fig. 6(b) shows the resulting locations of the triangle mesh vertices on the complex plane. The vertices coinciding with the hemispherical cut path are shown in red. The north pole point is mapped to the point at infinity of the extended complex plane.

Points on the complex plane can be mapped onto the unit sphere surface using the inverse stereographic projection, given by

where w²=u²+v², and u and v represent the real and imaginary coordinates of a point in the complex plane. This transformation yields points on the unit sphere (with the north pole removed). A simple way to picture this operation is to imagine a unit sphere sitting on the complex plane with its south pole on the origin and its north pole pointing straight up, out of the plane, as illustrated by Fig. 7. All points are mapped to the sphere by connecting the point itself to the north pole of the sphere with a line segment, and the point where the line segment intersects the sphere surface is the mapped point. The north pole of the unit sphere is identified with the point at infinity of the extended complex plane, for completeness. The inverse stereographic projection is conformal, so the combined process yields a conformal map from the cortex to the sphere.

Although it is possible to directly map the points on the complex plane produced by Angenent’s method to the sphere using inverse stereographic projection, it is also possible to take advantage of certain remaining flexibilities in the steps leading up to the inverse stereographic projection and in the projection process itself to yield more standard maps with smaller metric distortion and consistent orientation. Furthermore, we can move the sphere to which we map in order to yield separate hemispherical maps of the cortical hemispheres, each of which can be conformal to a partial sphere, yet having even smaller metric distortion than the complete spherical map. These steps are important in the production of a standard map that puts sulci and gyri from different brains into similar positions on the spherical or hemispherical maps. We will refer to either the spherical map or the hemispherical map as the final map when such generality is permitted in the following discussion.

Our spherical map is designed to map each hemisphere of the brain to one half of the sphere such that the hemispherical cut maps to a great circle, as shown in Fig. 6(c). It is also designed to be consistently oriented and appropriately scaled. By “consistently oriented”, we mean that the anterior part of the brain corresponds to the upper (positive y) part of the sphere and the posterior part of the brain corresponds to the lower (negative y) part of the sphere. Also, the left hemisphere should correspond to the left (negative x) part of the sphere and the right hemisphere should correspond to the right (positive x) part of the sphere. The map is also designed to distribute the points on the sphere such that the area distortion between the original cortical surface and the final map is minimized.

Two steps are required in order to map the hemispherical cut path to a great circle in the center of the unit sphere. First, we must define a south pole point on the brain along the hemispherical cut path, which will correspond to the south pole of the sphere. The south pole point is defined as the point on the hemispherical cut path with the maximum geodesic distance (cf. Rettmann et al., 2002) from the north pole point. All points on the complex plane are then translated so that the south pole point maps to the complex origin. Second, the points on the complex plane must be rotated such that the points to which the hemispherical cut path maps are oriented along the imaginary axis and anterior corresponds to the positive imaginary axis. This step guarantees that the anterior/ posterior orientation is correct. If necessary, in order to obtain the appropriate left/right orientation, the points on the complex plane are then reflected left-to-right to make the right hemisphere map to the half-plane corresponding to the positive real numbers. The distribution of the points on the sphere can be modified by uniformly scaling the magnitude of the points on the complex plane, leading to different stereographic projections. The one with the minimum area distortion with respect to the original surface (which is described below) will be the final map on the spherical surface.

Imagine rotating the unit sphere 90° about the y-axis so that the hemispherical cut path will be mapped to the equator of the unit sphere instead of a great circle through its equator, as illustrated in Figs. 8(a) and (b). This rotation facilitates a unified presentation of both spherical mapping and hemispherical mapping, in which the points on the hemispherical cut path fall on a horizontal circle within the sphere that is not necessarily the equator, as shown in Fig. 8(c). Defining θ as the angle off the z-axis, then this circle is defined by the angle any point on the circle makes with the z-axis, and the center of the circle is defined by the point (0,0,cosθ). For the spherical map, θ=90°, while in a hemispherical map, 0<θ<90°. With reference to Fig. 7, it is not hard to see that the radius of the hemispherical cut path in the complex plane should be r = 2sinθ/(1−cosθ) (centered at the origin), as shown in Figs. 8(e) and (f).

In a hemispherical map, the circumference of the circle to which the hemispherical cut maps is smaller than that on a spherical map. Also, one hemisphere of the cortex occupies a larger portion of the sphere than the other. It is intuitively clear that cortical distances will be better represented on a hemispherical map of the brain hemisphere that maps to the large portion of the sphere. It is logical that the area distortion would be smaller on this map as well. Accordingly, this larger map yields a hemispherical map for a given brain hemisphere, and the second hemispherical map can be computed by rotating the sphere 180° around the y-axis.

These steps, except the left-to-right reflection, can be undertaken using a single transformation in the complex plane called the Möbius transformation (Needham, 2000). A Möbius transformation is a linear fractional function M given by

where {a,b,c,d} are the coefficients of the Möbius transformation, ad−bc ≠0, {a, b, c, d} ∈ , and {z, z′} ∈ , where is defined as the complex plane. The transformation is analytic on \ {−d/c} and is one-to-one and conformal there. The transformation can be extended to the entire extended complex plane (Riemann sphere) ∪ {∞} by defining

Every Möbius transformation defines a one-to-one conformal mapping of a Riemann sphere to itself, and the composition of a conformal flattening with a Möbius transformation is still conformal. A Möbius transformation maps extended circles onto extended circles which is determined by three points and ∞ can be one of these points. Accordingly, it is completely determined by its values at three distinct points. The required translation, rotation and scaling can be used to select these three points, and the points they will be mapped to determine the a, b, c and d coefficients of the Möbius transformation we are interested in.

Consider the following particular Möbius transformation:

where {r} ∈ ℝ and {ρ, s} ∈ are the coefficients of the Möbius transformation, where ℝ is the set of real numbers. Eq. (6) performs a translation, rotation and scaling using the term ρ(z−s) (transforming Fig. 6(b) into Fig. 8(d), for example). Next, the unit sphere is rotated and scaled using the mapping r(z+1)/(z−1) (i.e., mapping the vertical line (on the imaginary axis) in Fig. 8(d) to a circle of radius r on the complex plane in Fig. 8(f)). With reference to Figs. 6 and 8, the relationship between the coefficients of the proposed Möbius transformation in Eq. (6) and the translation, rotation and scaling steps explained above are as follows: the coefficient s is determined by the required translation, the phase of (complex number) ρ is determined by the rotation, and the modulus of ρ and r are determined by the scaling required to minimize the area distortion by properly distributing the points on the unit sphere surface.

We defined a measure of area distortion of the final map to be

where T is the number of triangles in the mesh, B_i is the area of the ith triangle on the cortical surface, B is the total surface area on the cortical surface, S_i(|ρ|) is the area of the ith triangle on the sphere, which depends on the modulus of the scale factor ρ, and S(r) is the total surface area of the final map as a function of r,

The measure of area distortion D(r,|ρ|) is minimized within the group of Möbius transformations,

for a given s ∈ , using the conjugate gradient algorithm. r=2 yields a spherical map, where optimization on both hemispheres is carried out simultaneously. For hemispherical maps, separate optimization is carried out on the left and right cortical hemispheres. Based on our analysis of the relative area distortion between the final maps and the original cortical surfaces for fifteen human brains reported in (Tosun and Prince, 2001), the optimum conical angle of the spherical cap removed out of a total of 360° for both the left and right cortical hemispheres is 64° – i.e., θ=32°. The final maps on the spherical surface with segmented sulci (Rettmann et al., 2002) are shown in Fig. 9. Area distortion is substantially smaller and more uniform over the hemispherical maps as compared to the spherical map (Tosun and Prince, 2001).

3.3. Surface registration

Conformal mapping can be applied at any stage to map a surface to a sphere or hemisphere. Since triangle regularity of the meshes is quite good after surface inflation, the conformal mapping is less sensitive to its initial conditions. But it is still necessary to pick a “north pole” for the conformal mapping of a surface so that it corresponds to a common position across all brains – these north poles will map to the same point on the sphere.

To choose a good common north pole, we have selected one brain out of our database of 35 brains to be an “atlas brain”. All other partially flattened surfaces are registered to the partially flattened surface of our atlas brain using a modified iterated closest point (ICP) (Besl and McKay, 1992) algorithm. The north pole of the atlas brain, selected at the top of the hemispherical cut path is then projected to the subject brain and used as the subject’s north pole for conformal mapping. The surface registration provides not only a north pole point with a common position across all brains, but also standardizes the way the conformal coordinates are set up on the surfaces (see (Angenent et al., 1999)). Furthermore, registering the partially flattened surfaces of the atlas and the subject brains instead of their original surfaces permits the construction of the second spherical mapping method outlined in Fig. 3, and makes the method more consistent and stable than using the original surfaces which are highly variable and detailed compared to the partially flattened surfaces.

We use a modification of ICP to align a given partially flattened surface to the atlas’ partially flattened surface. ICP is an iterative method for registering two point sets. Given two surfaces V and W represented by the point sets ${\{v_i\}}_{i=1}^N$ in V and ${\{w_j\}}_{j=1}^M$ in W – mesh vertices connected by triangular surface elements – the ICP algorithm computes the rigid body transformation that aligns the point sets.

Ordinarily, ICP ignores the fact that V and W are surfaces when creating correspondences. Instead, it finds a mapping (i) that gives the closest vertex ∈ W corresponding to the vertex v_i∈V, forming the correspondence pair (v_i, ). At each iteration, ICP first finds the closest point correspondences between the point sets, and then finds the rigid body transformation (·) of points in W that minimizes the average distance of the corresponding points, given by

where ||·||₂ gives the length of a vector. These steps are repeated until convergence.

Our modification to ICP acknowledges the fact that V and W are triangulated surfaces rather than simple point sets. This means that the closest point to v_i∈V on W may not be a vertex, but may instead be a point on a triangular element defined by the vertices in W, as illustrated in Fig. 10. In our modification to ICP, correspondence pairs are identified by projection operators (·) (onto V) and (·) (onto W). In particular, we identify the closest point to v_i∈V on the surface W by the projection (v_i), and similarly, the closest point to w_j∈W on the surface V is given by (w_j). This identification yields the correspondence pairs (v_i, (v_i)), i=1,…,N and (w_j, (w_j)), j=1,…,M. In order to reduce possible directional bias, we define the modified ICP objective function in a bilateral fashion, as follows:

The fact that W is being moved under a rigid body transformation is reflected in this expression by the use of the rigid body transformation operator (·).

To provide quantitative comparison on the matching quality of ICP and the modified ICP algorithms, we defined two distance measures from one surface to the other. The first distance function “W-distance” is defined on the mesh nodes ${\{w_j\}}_{j=1}^M$ of the surface W, as the distance between correspondence pairs (w_j, (w_j)), j=1,…,M. In an analogous fashion, the second distance function “V-distance” is defined on the mesh nodes ${\{v_i\}}_{i=1}^N$ of the surface V, as the distance between correspondence pairs (v_j, (v_i)), i=1,…,N. W-distance and V-distance measures are illustrated in Fig. 10.

In order to assess the matching quality of the surface registration algorithms, each of the 34 subjects’ surfaces is registered to the surface of our atlas brain twice, once using the ICP algorithm and once using the modified ICP algorithm. The W-distances from the registered surfaces to the atlas surface and V-distances from the atlas surface to the registered surfaces are measured. For each subject, the modified ICP yields a combined distance (incorporating both W-distance and V-distance) with smaller mean and standard deviation over the mesh vertices, compared with those of the ICP algorithm. In Table 1, the range and the average of the means and standard deviations over 34 subjects are reported. Smaller mean values and significant decrease in the standard deviation support the claim that our modifications to the ICP algorithm improve the matching quality of the rigid body transformation.

3.4. Spherical mapping

The preceding alignment procedure automatically yields a common north pole and conformal coordinates so that a given partially flattened surface can generate its own conformal map to the sphere or hemisphere using the method described in Section 3.2. It is guaranteed that the orientation of the mapped surface is correct and that its area distortion relative to the original surface is minimized (within the group of Möbius transformations).

The first method we describe (see the block diagram shown in Fig. 3) starts with the iterative surface inflation stopped at the stage of a partially flattened map. Then, the partially flattened surface is aligned with the atlas’ partially flattened surface and the atlas’ north pole point is projected to the subject’s partially flattened surface. The maps on the spherical surface are then generated using the conformal mapping algorithm. There is no expectation, of course, that the total map will be conformal or isometric to the original cortical surface.

The second method (again, see the block diagram in Fig. 3) is readily developed using the single optimized conformal map of the atlas brain’s partially flattened surface. Since each subject’s partially flattened surface is registered to the atlas, and all surfaces have very similar overall shape due to the mean curvature stopping criterion, it is possible to map all subject’s partially flattened surface vertices to the atlas’ partially flattened surface. Then, since the atlas’ conformal map is a continuous bijection defined on the entire surface of its triangle mesh, all these points are readily mapped to the sphere (or hemisphere) using the atlas map, which is a pseudo-conformal map. In this way, to the extent that regions – e.g., sulci and gyri – are well-aligned on the partially flattened surfaces, they are guaranteed to be well-aligned on the spherical maps. This method proves to be quite reliable compared to the first method, as we demonstrate in the next section.

4.1. Probabilistic maps

We randomly picked 30 brain surfaces from our database of 35 surfaces. Four sulcal regions, defined as the buried cortex surrounding the sulcal spaces, were automatically segmented (Rettmann et al., 2002) and manually labelled on all 35 brains. These sulcal regions are the central sulcus (cs), superior frontal (sf), cingulate sulcus (cn) and parieto-occipital sulcus (po) on both the left and right cortical hemispheres. The reason we select these sulci is that they could be consistently labelled across three scans of the same subject (Rettmann, 2003). These sulci are colored on the left hemisphere of one cortical surface in Fig. 11.

Let HS_r denote the reference hemisphere to which all brains are mapped. (This hemisphere is parameterized by the particular conformal map of the atlas brain, which is just one of the 30 brains picked randomly.) To see where the four sulci map to on HS_r, we generated four probabilistic maps I_i:HS_r→[0,1] for i=cs, sf, cn and po. For each sulcal region, I_i is initialized to zero, and for each vertex of HS_r, the value of the probabilistic map is incremented by 1 for each training brain, if the corresponding point on the hemispherical map of that brain has the same sulcal label. Then I_i is divided by 30 to normalize its values to the range [0,1]. This gives a probabilistic map comprising an estimated conditional probability of the location of a sulcal region for each of four sulcal regions.

Two probabilistic maps were generated, one for the map I⁽¹⁾ created using individual conformal maps (first method), and one for the map I⁽²⁾ created by using only the atlas conformal map (second method). Fig. 12 shows the probabilistic maps $I_i^{(2)}$ for i=cs, sf, cn, po for both the left and right cortical hemispheres. The maximum values of the probabilistic maps $I_i^{(1)}$ and $I_i^{(2)}$ are given in Table 2. The table reveals a lack of good sulcal alignment for Method 1, the individual conformal maps method.

4.2. Automatic sulcal labelling

Here, we examine the performance of a simple automated sulcal labelling technique using the probabilistic maps derived in the previous section. Let I_i be one of the probabilistic maps derived in the previous section using the 30 training set surfaces. Now consider one the five remaining test surfaces, each of which is mapped to a hemisphere using either the first or second method. For each segmented sulcal region on the test surface, the maximum of each probabilistic map is calculated. If the maximum is larger than 0.5, then the label of that probabilistic map is assigned to the entire sulcal region.

Fig. 13(a) shows a mapped sulcal segmentation. The correct manual labels for the four sulcal regions under consideration are given (for the left hemisphere) in Fig. 13(b). Fig. 13(c) shows the result of the automated labelling technique.

To evaluate the effectiveness of this labelling technique, and hence the spherical mapping technique, two types of labelling errors are defined: probability of detection $P_i^{\mathrm{D}}$ and probability of false alarm $P_i^{\mathrm{F}}$. Probability of detection is defined as the percentage area of the manually labelled region that is correctly labelled by our automatic labelling method. Probability of false alarm is defined as the percentage area of the automatically labelled region that is not manually labelled. The probabilities of detection and false alarm averaged for the five test brains are reported in Table 3. Higher $P_i^{\mathrm{D}}$ values and lower $P_i^{\mathrm{F}}$ show that our second method yields a better standardized coordinate system. The probability of detection in the range between 0.83–1.00 is quite good in comparison to other methods in the literature (cf. Behnke et al., 2003).

4.3. Quality of isometry

In this work, our aim is to develop a standard, automatic mapping technique that makes comparison across multiple brains possible. Although moderate distortion of distance can be achieved in spherical mapping techniques, it is a well known fact (Carmo, 1976) that a cortical surface and a spherical surface are not developable – i.e., there is no isomorphism between a cortical surface and a sphere since a cortical surface has spatially varying Gaussian curvature unlike a sphere’s constant Gaussian curvature.

To quantify the degree of isometry between a cortical surface and its final mapped surface, a random reference vertex X and a set of N random vertices ${\{Y_i\}}_{i=1}^N$ are picked on the given cortical surface. Then the minimal geodesic distance from X to each Y_i on the original surface and on the final map is calculated by the fast marching algorithm on the triangle mesh (cf. Rettmann et al., 2002), and denoted by d_o(i) and d_f(i), respectively. To avoid possible errors due to distance normalization required for comparison, the pairwise ratio of geodesic distances is used instead of absolute geodesic distance. Let the pairwise ratio of geodesic distances ${\{d_o(i)\}}_{i=1}^N$ be organized as a lower triangular matrix T, T_ij=d_o(i)/ d_o(j), and the pairwise ratio of geodesic distances ${\{d_f(i)\}}_{i=1}^N$ be organized as a lower triangular matrix T^*, $T_{ij}^{\ast }=d_f(i)/d_f(j)$. A similarity measure between T and T^* is defined by E:

similar to the “stress” measure proposed by (Sammon, 1969; Schwartz et al., 1989) to quantify the error in the “multidimensional” scaling problem. To eliminate the dependency of measure E on the reference point selection, for each of 35 subjects, E is measured 100 times using 100 different reference points X, and point sets ${\{Y_i\}}_{i=1}^N$ of size N=5000. The analysis of variance (ANOVA) technique (Scheffe, 1959; Crow et al., 1990) is used to decide whether different reference point selections have a significant effect on measure E for a fixed value of N=5000, and to estimate these effects.

Statistics over 35 subjects are reported in Table 4. The reference point selection factor failed to reach significance (Pr(>F) >0.1) for any combination of original surfaces and final maps analyzed in this experiment. The overall mean stress value on the conformal mapping and pseudo-conformal is between 0.044 and 0.070 which are comparable with the errors reported in (Sammon, 1969). Note that a lower value of E indicates a better level of isometry. As expected, a better level of isometry is achieved in the hemispherical maps as compared to the spherical maps. For partially flattened surfaces, the mean stress value is 0.013 which is comparable to the stress value about 0.010 achieved by the cortical surface flattening method proposed in (Schwartz et al., 1989). The reason for the better isometry level of the partially flattened surfaces is due to the geometric features – i.e., deep sulcal regions – preserved on these maps that the spherical maps lack.