Hierarchical spherical deformation for cortical surface registration

Abstract

We present hierarchical spherical deformation for a group-wise shape correspondence to address template selection bias and to minimize registration distortion. In this work, we aim at a continuous and smooth deformation field to guide accurate cortical surface registration. In conventional spherical registration methods, a global rigid alignment and local deformation are independently performed. Motivated by the composition of precession and intrinsic rotation, we simultaneously optimize global rigid rotation and non-rigid local deformation by utilizing spherical harmonics interpolation of local composite rotations in a single framework. To this end, we indirectly encode local displacements by such local composite rotations as functions of spherical locations. Furthermore, we introduce an additional regularization term to the spherical deformation, which maximizes its rigidity while reducing registration distortion. To improve surface registration performance, we employ the second order approximation of the energy function that enables fast convergence of the optimization. In the experiments, we validate our method on healthy normal subjects with manual cortical surface parcellation in registration accuracy and distortion. We show an improved shape correspondence with high accuracy in cortical surface parcellation and significantly low registration distortion in surface area and edge length. In addition to validation, we discuss parameter tuning, optimization, and implementation design with potential acceleration.

Keywords: Cortical surface registration; Shape correspondence; Spherical deformation; Spherical harmonics.

Figure 1:

An example of angular interpolation failures. The displacements by a counter clockwise rotation about the fixed axis (blue) are represented by elevation angles passing through the pole (green). Such angles have different signs before and after the pole (red and purple). The resulting interpolation thus yields rotation singularity at the pole, which is incapable of encoding the rotation completely.

Figure 2:

Precession and intrinsic rotation: (a) initial setting of two frames, (b) z-axis alignment after precession, and (c) the final alignment after intrinsic rotation. Any rigid rotation can be implemented by precession and intrinsic rotation. The resulting composite rotation does not rely on a particular spherical coordinate system.

Figure 3:

A schematic illustration of the proposed rotation by the axis-angle representation. For the rotation of a given location, (precession) the rotation axis z (red) is rotated to $\widehat{\mathbf{\text{z}}}$ (blue) by ω^⊥ about z^⊥, followed by (intrinsic rotation) a rotation about $\widehat{\mathbf{\text{z}}}$ by ω (green). The exponential map (purple) at z is employed to encode local geodesics (orange). Finally, the rotation axis $\widehat{\mathbf{\text{z}}}$ and its associated rotation angle ω smoothly vary on the unit sphere as functions of spherical locations. A half sphere is used for better visualization.

Figure 4:

The average hCurv feature maps at intermediate spherical harmonics degrees l (optimization with single resolution based on only hCurv). Each hemisphere shows the average hCurv feature after independent optimization at each individual degree. The cortical folding patterns become sharper, and the finest patterns are achieved after spherical harmonics coefficients are optimized together.

Figure 5:

The average hCurv feature maps on the 30 subjects after pair-wise registration to a fixed template. Overall, these methods achieve similar hCurv patterns. The proposed method including a non-optimal rigid alignment provides a sharper representation close to the template.

Figure 6:

An example of pair-wise registration from a single subject to a fixed template. These methods begin with almost the same rigid alignment before the local deformation and produce similar hCurv patterns with slight difference. However, FreeSurfer and Spherical Demons only update the local deformation during the optimization. This results in relatively large deformation in several regions, whereas the proposed method updates both rigid and non-rigid deformation to reduce locally focused deformation. Also, the optimal rigid alignment in our method provides improved feature alignments and registration distortion compared to ours with a non-optimal rigid alignment.

Figure 7:

The average hCurv feature maps on the 30 subjects. The three methods achieve similar hCurv patterns, while FreeSurfer shows little more blurred patterns than other methods. These methods provide much more improved average patterns than the initial average.

Figure 8:

hCurv variance at different degree l. The feature variance decreases as l increases. Both hemispheres have similar variance at each degree. The proposed method has smaller variance then FreeSurfer and Spherical Demons at l ≥ 9 and l ≥ 19, respectively. It is noteworthy that smaller variance does not necessarily indicate better surface registration performance. In our experiments, the proposed method works well at l = 15 in terms of cortical parcellation and registration distortion.

Figure 9:

Dice coefficient of 49 regions on the left and right hemispheres. One-sided t-tests reveal regions with statistical significance after the FDR correction (q = 0.05). Several regions are significantly improved, while no region becomes worse after the FDR correction. The color in the labels indicates the improved regions compared to FreeSurfer (blue), Spherical Demons (red), and both methods (green). In comparison with FreeSurfer, our method has 21 and 12 improved regions for the left and right hemispheres, and with Spherical Demons, our method has 9 and 7 improved regions for the left and right hemispheres (see Fig. 10 for the improved regions with the adjusted p-values).

Figure 10:

Negative log of the adjusted p-values on cortical regions with significantly improved Dice coefficients after the FDR correction (q = 0.05). The average surface is divided by the mode map of 49 regions. Total 33 and 16 out of 98 regions are significantly improved compared to FreeSurfer and Spherical Demons, respectively. The color indicates negative log of the adjusted p-values in the improved regions.

Figure 11:

Area distortion (whiskers with maximum 1.5 interquartile range). The proposed method has a less skewed distribution to the right (shorter tail). This implies a fewer number of regions with large area distortion than FreeSurfer and Spherical Demons on both hemispheres.

Figure 12:

Edge distortion (whiskers with maximum 1.5 interquartile range). The proposed method has a less skewed distribution to the right (shorter tail). This implies a fewer number of regions with large edge distortion than FreeSurfer and Spherical Demons on both hemispheres.

Figure 13:

Area change of 49 regions on the left and right hemispheres. One-sided t-tests reveal regions with statistical significance after the FDR correction (q = 0.05). More than a one third of regions have significantly reduced area change, while no region becomes worse after the FDR correction. In comparison with FreeSurfer, our method has 18 and 21 improved regions for the left and right hemispheres, and with Spherical Demons, our method has 16 and 23 improved regions for the left and right hemispheres (see Fig. 14 for the improved regions with the adjusted p-values). Note that the maximum range is truncated at 60% for better visualization. The color in the labels indicates the improved regions compared to FreeSurfer (blue), Spherical Demons (red), and both methods (green).

Figure 14:

Negative log of the adjusted p-values on cortical regions with significantly reduced area change after the FDR correction (q = 0.05). The average surface is divided by the mode map of 49 regions. Total 39 out of 98 regions are significantly improved compared to FreeSurfer and Spherical Demons. The color indicates negative log of the adjusted p-values in the improved regions.

Figure 15:

An example area change in Ent. (top) The three methods yield similar cortical folding patterns after co-registration (see Fig. 7 for the average hCurv maps). The mode regions of Ent are highlighted brightly. (middle) The Dice coefficients are comparable, which implies that these methods achieve comparable performance in surface alignment; the hCurv maps are also well aligned with the averages. (bottom) Even with comparable registration performance, the surface area (triangle size) is less distorted in the proposed method than FreeSurfer and Spherical Demons. It is noteworthy that the mode region of Ent in our group-wise framework is little larger than others because the distortion can be better minimized (i.e., better area preservation) in this way, while maintaining comparable registration accuracy (see Fig. 9 for Ent).