Establishment of anatomical and functional correspondence is a crucial step toward gaining insight into biological processes. Neuroscience research efforts, such as characterizing brain morphology, require accurate and robust methods for producing such mappings. The extensive literature detailing methodology is evidence of the rich history of algorithmic development which continues contemporaneously. We highlight several key historical contributions which are particularly relevant to the work presented.

Free-form deformation (FFD) image registration, characterized by regularization based on the B-spline basis functions, has several advantages including algorithmic simplicity, good performance, and guaranteed parametric continuity. Current research was preceded by related work for geometric modeling (Sederberg and Parry, 1986) and originated with such important contributions as Szeliski and Coughlan (1997); Thévenaz et al. (1998), and Rueckert et al. (1999). Continued development within this early spline-based paradigm produced additional innovations such as integrated similarity metrics (e.g., Mattes et al., 2003), additional transformation constraints (e.g., Rohlfing et al., 2003), and notable open source implementations (e.g., Ibanez et al., 2005; Klein et al., 2010b; Modat et al., 2010; Shackleford et al., 2010).

Parallel to this branch of algorithmic progress are the informally denoted “dense transforms” perhaps best exemplified by Thirion's seminal contribution (Thirion, 1998). Relationships with earlier elastic (Bajcsy and Kovacic, 1989; Gee et al., 1993) and fluid (Christensen et al., 1996) registration methods are detailed in the works of Bro-Nielsen and Gramkow (1996) and Pennec et al. (1999) who observe that smoothing via Gaussian convolution, a defining characteristic of Demons, of the update or total displacement field is a greedy approximation for solving the partial differential equations governing the physics of an elastic or fluid deformation, respectively. However, the use of such approximations entails that physical properties, such as topological regularity, are no longer guaranteed.

It is interesting to note that within this context, traditional FFD algorithms can be viewed as a type of fluid-like Demons approach where, rather than projecting the update field to the space of regularized fields using Gaussian convolution, gradient fields are projected to a smooth space characterized by the B-spline basis functions. This analogy was hinted at in our earlier work (Tustison et al., 2009) where we showed that fitting the update field to a B-spline object using a fast approximation routine (Tustison and Gee, 2006) is equivalent to a preconditioning of the standard gradient used in gradient descent-based FFD optimization. This preconditioning is used to mitigate the hemstitching effect induced by the ill-conditioned nature of the traditional gradient-based FFD formulation¹. We denoted this new FFD variant as directly manipulated free-form deformation (DMFFD) and, as part of the ITKv4 refactoring efforts, has been implemented for use with the new registration framework² which permits both B-spline smoothing on the update (“viscous”) and total (“elastic”) displacement fields at each iteration (cf analogous Gaussian, i.e., Demons, implementation³).

Continuing from the work of Christensen et al. (1996) and subsequent exploration into the mathematical formalisms of diffeomorphisms (e.g., Dupuis and Grenander, 1998), the well-known Large Deformation Diffeomorphic Metric Mapping (LDDMM) algorithm was proposed in Beg et al. (2005). In contrast to the mapping produced by Christensen et al. (1996), LDDMM yields the geodesic solution in the space of diffeomorphisms between two images. Since its introduction, LDDMM has inspired much innovation in the image registration literature. Applying the log-Euclidean framework of Arsigny et al. (2006), DARTEL (Diffeomorphic Anatomic Registration using Exponential Lie algebra) uses a constant velocity field parameterization to provide a fast, diffeomorphic alternative (Ashburner, 2007). Additionally, symmetrical considerations in the velocity field parameterization are discussed in Avants et al. (2008) in the context of a cross correlation similarity metric. By explicitly symmetrizing the LDDMM formulation, this Symmetric Normalization (SyN) approach minimizes the bias of the resulting transformation when selecting the “fixed” and “moving” images. A greedy version of this algorithm has proven successful in neuroimaging (Klein et al., 2009) and pulmonary (Murphy et al., 2011) applications as well as in multi-atlas label fusion (Wang et al., 2013).

Although many extensions of LDDMM rely on some form of Gaussian convolution for regularization (e.g., Risser et al., 2011), there has been significant interest in constraining FFD approaches to the space of diffeomorphisms. An early attempt reported in Rueckert et al. (2006) enforced diffeomorphic transforms by concatenating multiple FFD transforms, each of which is constrained to describe a one-to-one mapping. Modat et al. (2011) incorporated the log-Euclidean framework for enforcing diffeomorphic transformations and ensuring invertibility. Similarly, the work of De Craene et al. (2011) provided a full LDDMM-style algorithm based on B-splines called temporal FFD in which the time-varying velocity field is modeled using a 4-D B-spline object (3-D + time). Numerical Eulerian integration of the mapping propagated within the velocity field yields the transform between parameterized time points.

As alluded to earlier, B-spline approximation can also be used for regularizing time-varying vector fields in an analogous fashion as Gaussian convolution. In this vein, and similar to De Craene et al. (2011), we reported in Tustison and Avants (2012a,b) the use of an n-D + time B-spline object to represent the characteristic velocity fields. However, we use the DMFFD formulation to improve the solution convergence. This also facilitates modeling temporal periodicity and the enforcement of stationary boundaries. Both this work and our earlier work (Tustison et al., 2009) demonstrate that the DMFFD framework is potentially applicable to the entire gamut of diffeomorphic registration algorithms and provides alternative smoothing possibilities with different continuity properties (e.g., C² vs. C³).

The two regularization approaches (i.e., Gaussian convolution vs. DMFFD), however, produce characteristically different solutions. In addition to smoothing kernel differences, Gaussian convolution tends to “flatten” the signal in contrast to an approximation or fitting of the signal provided by the DMFFD approach. Also, whereas Gaussian convolution operates entirely within discretized space, the B-spline approximation routine constructs a continuous object prior to any voxelwise reconstruction of the sampled fields. Interestingly, a similar comparison was made with respect to Gaussian derivative estimation (Bouma et al., 2007). Although typically estimated using truncated, discrete Gaussian convolution, an alternative based on B-spline approximation demonstrated superior performance with similar computational cost.

Three popular diffeomorphic algorithms and their DMFFD analogs (LDDMM, DARTEL, and SyN) were implemented by the authors as part of the recent refactoring of the open source Insight Toolkit (ITK) although related work had been previously implemented within the popular Advanced Normalization Tools (ANTs)⁴. Given the popularity and excellent performance of the greedy variant of SyN, our evaluation focus in this work is its B-spline analog which we denote as “B-spline SyN” or “DMFFD SyN.” Evaluations of the respective algorithmic instantiations are performed using the antsRegistration program found in the ANTs repository (also originally developed by the authors). This permits a direct algorithmic comparison as potential sources for implementation bias have been reduced (Tustison et al., 2013). Additionally, in the spirit of open science, all text, figures, and scripts to reproduce the results contained in this work are publicly available online⁵.

As mentioned previously, a template was constructed for each data set (cf Figure 2) from all cohort images. Subsequently, each image was registered to its corresponding template using either SyN or B-spline SyN as described previously (prior linear registration stages were identical between the two algorithms). As a brute-force parameter exploration is not a part of this work, we rely on previously reported research (Klein et al., 2009; Avants et al., 2011) and our own experience as authors/developers of the algorithm/software to select parameters which demonstrate robust performance across data sets. For both algorithms, the gradient step was 0.1 for each of the four multi-resolution levels with shrink factors of {6, 4, 2, 1} and Gaussian smoothing for each of those levels being (0, {9, 4, 1, 0}) in terms of voxels. The number of iterations per level were {100, 100, 70, 20} with a convergence threshold of 10⁻⁹ and window size of 15 iterations²¹.

All processing was performed using the linux cluster at the University of Virginia²² using the PBS Pro queuing system for managing resources. The perl scripts used to create the jobs for the cluster are included in the github account associated with this evaluation. For the data sets used in this study, times for B-spline SyN were approximately 15–40% greater than Gaussian-based SyN using single-threading and a dense metric gradient sampling. Timing data for specific data sets are given in Table 2 Timing includes both rigid and affine transform optimizations.

The only difference between the two registration settings consists of the Gaussian and B-spline parameters governing the update field smoothing, S_v. In our experience, smoothing of the total field did not improve the results, at least for these data (which conforms with our experience with other data), so the total field smoothing, S_ϕ is 0 for both registration approaches. Specifically, the chosen parameters for the SyN algorithm were: S_ϕ = (0, 0) and S_v = (0, 3) in voxel terms²³. Although our experience with B-spline SyN is much more limited, we were able to choose comparable parameters based on a knot spacing for the update field of 26 mm at the base level which is reduced by a factor of two for each succeeding multiresolution level. This yields a final knot spacing of 3.25 mm²⁴. For comparison, after selecting these smoothing parameters we discovered in the supplementary material of Klein et al. (2009) the similarity to the gradient smoothing parameter for the IRTK FFD algorithm which also used four multi-resolution levels with an initial knot spacing of 20 mm per dimension for a final knot spacing of 2.5 mm.

Quality of overlap using the Dice similarity metric was determined from the transformed labels using the open source ITK implementation described in Tustison and Gee (2009). Both the original labels and MALF labels were warped to the fixed image for comparison using nearest neighbor interpolation. A joint Dice metric value was calculated from the combined labels for each cohort for each of the two SyN methods. These values are rendered in notched box plot format in Figure 5. Non-overlapping notches indicate approximately statistically significantly different median values at the 95% confidence level (McGill et al., 1978). For all data sets, the B-spline SyN variant showed a small but statistically significant improvement in overall Dice values. In order to provide a more complete picture of performance differences, we also accounted for label volumetric considerations (Rohlfing, 2012). In Figure 5 we plotted the Dice value difference between each SyN variant (B-spline SyN—SyN) for each label of each intra-subject registration pair within a data set vs. the volume of the label in the fixed image. Values above and below the dashed line indicate better regional performance for B-spline SyN and SyN, respectively.

To further characterize the deformable transform differences, we calculated the log of the Jacobian determinant of the transformations from each subject to the template and tabulated statistical information within the brain region only. A noticeable difference between the two algorithms was the respective range of values in log Jacobians. We plotted the (95th%—5th%) for each algorithm across all data sets in Figure 6. It is apparent that B-spline SyN results exhibit a much greater range of deformation. Qualitative differences are shown in Figure 7. In order to ensure randomness to minimize presentation bias in illustrating qualitative results (and given the relative poor performance in humans as random number generators Wagenaar, 1972), we used R to generate uniform random numbers for both subject and axial slice selection. We then used the log Jacobian images to locate regions of maximal difference between the SyN and B-spline SyN results. From Figure 7 it is quite apparent that the results are very similar which is to be expected considering the almost identical algorithmic make-up between the two approaches (e.g., similarity metric, implementation, linear transforms). However, there are subtle differences particularly in the cortex which help explain both the relative difference in Jacobian and Dice distribution.

B-spline SyN produced slightly greater Dice values than the original SyN. Although actual differences are relatively small, they are statistically significant. By implementing both algorithms in the same code base, we are not only able to eliminate non-regularization components of the registration but we are also able to eliminate implementation differences. Thus, performance disparity can be isolated to smoothing choice. However, even with this restricted focus there are various reasons for the evaluation outcome. These include the approximation-vs.-convolution distinction mentioned earlier for the two regularization approaches (which could also explain the reason why the range of log Jacobian values tend to be significantly higher for B-spline SyN). Also, the fact that the regularization for B-spline SyN is theoretically continuous whereas the truncated Gaussian convolution is only a discrete approximation could be a potential factor.

Additional observations of interest concern the differences in results between the MALF and original labelings for all data sets. Not only was there an overall increase in performance for both algorithms with the MALF labels for all cohorts, but there was also an increase in performance disparity relative to the variance in the resulting Dice values. A possible explanation for this, and one that seems quite plausible, is that the MALF labelings are derived from registering a cohort to the target subject and then performing a consensus labeling (Wang et al., 2013). Both these steps are heavily reliant on image intensity information (in fact, both use a form of correlation as the similarity measure for optimization). Since the MALF “correction” tends to group labeled regions according to the same metric used for establishing anatomical correspondences, alignment of these labeled regions seems much more likely which would result in higher Dice metrics. A related effect for labeled data in general is being currently investigated by the authors. From the label volume vs. Dice difference plots, there is no immediately discernible pattern of performance variation with region size. However, a possible confound is that region definitions vary between cohorts. Although we did not look at region vs. performance difference variation within cohorts, such inquiry is certainly possible as we have made the resulting csv files available with the github repository associated with this work.

Relative to Klein's study (Klein et al., 2009), it should be noted that the SyN implementation used in that study is found in the ANTS program which is a precursor of the antsRegistration program described earlier. Although the theoretical aspects are the same, there are substantial implementation differences between the two programs. In addition, several parameters varied between the two studies which translated into a more aggressive metric and gradient step in addition to fewer levels and iterations²⁵. In this study we took a more conservative approach based on our continued development and experience resulting in parameters which have proven useful in our cortical thickness pipeline (encapsulated in the ANTs script antsCorticalThickness.sh) and our experience with the MICCAI 2013 SATA challenge data. We also ran our own internal experiments with Gaussian-based SyN (both ANTS and antsRegistration) using the same conservative parameters on the MAL35 data for which the latter demonstrated slightly improved performance over the former.

Despite the thorough evaluation with multiple data sets, we readily acknowledge the limitations of this study including a very focused application, i.e., healthy brains of a single modality, and absence of a thorough exploration of parameter selection and sensitivity. Although such work might be beneficial (e.g., by aiding other researchers in parameter selection), characterizing parameter permutations of potential interest would expand the current work far beyond its intended scope. However, in addition to this work, SyN has also been previously evaluated in Klein et al. (2009) and Avants et al. (2011) and, based on additional experience and application, the parameter set was modified each time but still yielded excellent performance providing evidence for flexibility in parameter selection. Outside of a range of parameters based on sound engineering principles, experience, and intimate knowledge of both the corresponding algorithms and software, determination of optimal generic parameters even for a specific application is difficult. In fact, the “No Free Lunch Theorem” Wolpert and Macready (1997) emphasizes the importance of prior knowledge in tuning optimization algorithms for a particular application.

One of the advantages that has not been explored in this work is the use of B-spline SyN for small deformation estimation problems such as in pulmonary or cardiac applications. Such problems typically require greater regularization which implies larger discrete kernels for Gaussian convolution. A related issue concerns applications involving severely anisotropic data where the continuous nature of the DMFFD approach might help over Gaussian convolution. Also, we emphasize that underlying the DMFFD approaches is a fitting routine for sparse and scattered data which offers added flexibility over smoothing using discrete convolution where the latter implies regularly placed data on a rectilinear grid for conventional implementations. This advantage could translate into faster running times if only select points are used to drive the registration or make possible more complex registration scenarios involving data arranged continuously within a finite domain (e.g., Tustison et al., 2011). Finally, the possibility of varying data confidence values, as introduced in Tustison and Gee (2006) with DMFFD-based routines would permit incorporation of spatial preferential weighting (i.e., additional prior knowledge) during optimization. Ongoing work will continue to explore these issues.

A significant amount of research has been devoted to image registration algorithmic development. Given their many salient characteristics particularly with respect to large deformation estimation constrained by topological continuity, diffeomorphic registration approaches have been a particular focus in the neuroimaging community. However, many groups continue to find success with non-diffeomorphic FFD methods (e.g., Rueckert et al., 1999; Klein et al., 2010b). Using our DMFFD framework, B-spline regularization is easily adapted into the diffeomorphic registration framework and performs well compared to analogous algorithms which we demonstrated in this work for the case of the widely-used SyN.