Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

The Insight ToolKit image registration framework

The Insight ToolKit image registration framework ORIGINAL RESEARCH ARTICLE published: 28 April 2014 NEUROINFORMATICS doi: 10.3389/fninf.2014.00044 1 2 1 1 1 Brian B. Avants *, Nicholas J. Tustison , Michael Stauffer , Gang Song , Baohua Wu and James C. Gee Penn Image Computing and Science Laboratory, Department of Radiology, University of Pennsylvania, Philadelphia, PA, USA Department of Radiology and Medical Imaging, University of Virginia, Charlottesville, VA, USA Edited by: Publicly available scientific resources help establish evaluation standards, provide a Hans J. Johnson, The University of platform for teaching and improve reproducibility. Version 4 of the Insight ToolKit ( ITK ) Iowa, USA seeks to establish new standards in publicly available image registration methodology. Reviewed by: 4 4 ITK makes several advances in comparison to previous versions of ITK. ITK supports Marcel Prastawa, University of both multivariate images and objective functions; it also unifies high-dimensional Utah, USA Andrey Fedorov, Brigham and (deformation field) and low-dimensional (affine) transformations with metrics that are Women’s Hospital, USA reusable across transform types and with composite transforms that allow arbitrary *Correspondence: series of geometric mappings to be chained together seamlessly. Metrics and optimizers Brian B. Avants, Department of take advantage of multi-core resources, when available. Furthermore, ITK reduces the Radiology, University of parameter optimization burden via principled heuristics that automatically set scaling Pennsylvania, 3600 Market St., Suite 370, Philadelphia, PA 19104, USA across disparate parameter types (rotations vs. translations). A related approach also e-mail: avants@grasp.cis.upenn.edu constrains steps sizes for gradient-based optimizers. The result is that tuning for different metrics and/or image pairs is rarely necessary allowing the researcher to more easily focus on design/comparison of registration strategies. In total, the ITK contribution is intended as a structure to support reproducible research practices, will provide a more extensive foundation against which to evaluate new work in image registration and also enable application level programmers a broad suite of tools on which to build. Finally, we contextualize this work with a reference registration evaluation study with application to pediatric brain labeling. Keywords: registration, MRI, brain, open-source, death 1. INTRODUCTION resources have become critical to setting performance standards in international challenges that evaluate “real world” registra- As image registration methods mature—and their capabilities tion scenarios (see, for instance, the SATA 2013 and BRATS 2013 become more widely recognized—the number of applications challenges at MICCAI in Nagoya, Japan). increase (Rueckert et al., 1999; van Dalen et al., 2004; Miller et al., 2005; Shelton et al., 2005; Chen et al., 2008; Baloch and One source of benchmark methodology is the Insight ToolKit (ITK) (Yoo et al., 2002; Ackerman and Yoo, 2003), which marked Davatzikos, 2009; Cheung and Krishnan, 2009; Peyrat et al., 2010; Fedorov et al., 2011; Kikinis and Pieper, 2011; Metz a significant contribution to medical image processing when it et al., 2011; Murphy et al., 2011). Consequently, image regis- first emerged at the turn of the millennium. Since that time, ITK tration transitioned from being a field of active research, and has become a standard-bearer for image processing algorithms and, in particular, for image registration methods. In a review few applied results, to a field where the main focus is trans- lational. Image registration is now used to derive quantita- of ITK user interests, image registration was cited as the most important contribution of ITK (personal communication with tive biomarkers from images (Jack et al., 2010), plays a major role in business models and clinical products (especially in Terry Yoo). Numerous papers use ITK algorithms as standard ref- erences for implementations of Demons registration and mutual radiation oncology) (Cheung and Krishnan, 2009), has led to numerous new findings in studies of brain and behavior (e.g., information-based affine or B-Spline registration (van Dalen Bearden et al., 2007) and is a critical component in applica- et al., 2004; Shelton et al., 2005; Floca and Dickhaus, 2007; Chen tions in pathology, microscopy, surgical planning, and more et al., 2008; Cheung and Krishnan, 2009). Multiple toolkits extend (Miller et al., 2005; Shelton et al., 2005; Floca and Dickhaus, ITK registration methods in unique ways. Elastix provides very 2007; Chen et al., 2008; Cheung and Krishnan, 2009; Peyrat fast and accurate B-Spline registration (Klein et al., 2010; Murphy et al., 2010; Kikinis and Pieper, 2011; Murphy et al., 2011). et al., 2011). The diffeomorphic demons is a fast/efficient approx- Despite the increasing relevance of image registration across imation to a diffeomorphic mapping (Vercauteren et al., 2009). application domains, there are relatively few reference algorithm ANTs provides both flexibility and high average performance (Avants et al., 2011). The BRAINSFit algorithm is integrated into implementations available to the community. Furthermore, these Slicer for user-guided registration (Kikinis and Pieper, 2011). Each of these toolkits has both strengths and weaknesses (Klein This work is supported by National Library of Medicine sponsored ARRA stimulus funding. et al., 2010; Murphy et al., 2011) and was enabled by an ITK core. Frontiers in Neuroinformatics www.frontiersin.org April 2014 | Volume 8 | Article 44 | 1 Avants et al. ITKv4 image registration The Insight ToolKit began a major refactoring effort in 2010. optimization technique (in the Optimizersv4 directory), the The refactoring aimed to both simplify and extend the tech- metric measuring the registration quality (the Metricsv4 niques available in version 3.x with methods and ideas from a directory), the images or other data objects that enter the met- new set of prior work (Christensen et al., 1996; Rueckert et al., ric and the parameters that are being optimized. The parame- 1999; Jenkinson and Smith, 2001; Miller et al., 2005; Peyrat et al., ters are usually defined by a geometric transformation but may 2010; Avants et al., 2011). To make this technology more accessi- point to other relevant objects. Any of ITK ’s transforma- ble, ITK unifies the dense registration framework (displacement tions may be optimized by the framework. New transformations, field, diffeomorphisms) with the low-dimensional (B-Spline, relative to ITK ,include the DisplacementField trans- affine, rigid) framework by introducing composite transforms, forms that are useful for engendering Demons or B-Spline reg- deformation field transforms, and specializations that allowed istration strategies. New VelocityField transforms are also these to be optimized efficiently. A sub-goal set for ITK was available. A typical application developer would employ all of to simplify parameter setting by adding helper methods that these components. A good starting point for new users who use well-known principles of image registration to automatically wish to see how these tools work together, in source code, is scale transform components and set optimization parameters. foundinthe tests. Forinstance,see thefiles itkTimeVaryin ITK transforms are also newly applicable to objects such as vec- gBSplineVelocityFieldImageRegistrationTest. tors and tensors and will take into account covariant geometry cxx for an example of a B-Spline diffeomorphism applica- if necessary. Finally, ITK reconfigures the registration frame- tion, itkSyNImageRegistrationTest.cxx to see SyN in work to maximize multi-threading resources where possible. The ITK and itkSimpleImageRegistrationTest2.cxx revised registration framework within ITK is more thoroughly for a more basic example. integrated across transform models, is thread-safe and provides Several usability goals spurred ITK development. We sum- broader functionalitythaninprior releases. marize these here. David Donoho once commented (in paraphrase) that aca- demic publications are merely “advertisements” for the real work 2.1.2. Image registration should be achievable in one step which is constituted by the “complete instruction set” that pro- This overarching goal is best illustrated by Registration duces the results reported in the publication (Buckheit and Methodsv4 in which a user may string together a series Donoho, 1995). The first part of the remainder of this docu- of registration tools to perform (for instance) a transla- ment will provide an “advertisement” for the ITK framework and tion registration, followed by an affine registration, followed 3 4 summarize its evolution from ITK to ITK . We then detail by a diffeomorphic mapping each of which might use a potential applications of this ITK framework in the context different image similarity metric. The different transforms of a general nomenclature. While this work is indeed incom- are accumulated in the new itkCompositeTransform plete, in the sense of Donoho, we refer to source code and data which chains transforms together as in Figure 2. Thus, this when relevant. Furthermore, section 3.1 shows a series of repro- framework provides unprecedented ability to perform com- ducible examples of ITK in action. Several areas relevant to plex and staged registration mappings. Furthermore, the neuroinformatics are highlighted in these examples: optimal tem- frameworks automated parameter scaling, itkRegistrat plate construction, “challenging” registration scenarios involving ionParameterScalesEstimator, vastly reduces the dif- brain mapping in the presence of lesions or resection, registra- ficulty of tuning parameters for different transform/metric tion when initialization priors are weak, asymmetry analyses, combinations. functional MRI, and non-traditional registration strategies are all highlighted. We also establish performance benchmarks for the 2.1.3. ITK Transforms should be unified current ITK registration, in comparison to a method developed Each ITK transform now has either global support (affine trans- for ITK , via a standard brain labeling task. Finally, we discuss form) or local (or compact) support (a displacement field trans- future developments in the framework. form). If any map in a composite transform has global support then the composite transform has global support. Both “fixed” 2. MATERIALS AND METHODS and “moving” images may have initial transforms. This allows one 2.1. OVERVIEW OF THE UNIFIED FRAMEWORK to reduce “registration bias” that may be induced by asymmetric interpolation (Yushkevich et al., 2010). The overall purpose of the registration refactoring for ITK was to simplify the user experience and to accelerate and improve per- 2.1.4. Registration mappings should be applicable to a number of formance. Here, we summarize how ITK works toward these goals. popular data types, including DTI Our revisions to the ITK transform hierarchy validated and extended the ITK transforms for thread safety and applicability 2.1.1. Core Software Components to not only vectors but also tensors. Reorientation steps necessary Figure 1 sketches the ITK architecture at a high level. for diffusion tensor mappings are now included in ITK . Registration applications are known as “registration methods” as they were in ITK . The methods, with source contained in ITK ’s RegistrationMethodsv4 directory, hold together 2 4 All ITK metrics are set to be minimized. For instance, the the different subcomponents that make a working instantia- itkMattesMutualInformationImageToImageMetricv4 returns tion of a registration strategy. These subcomponents include the negative mutual information. Frontiers in Neuroinformatics www.frontiersin.org April 2014 | Volume 8 | Article 44 | 2 Avants et al. ITKv4 image registration ITKv4 Registration Framework FIGURE 1 | A schematic overview of the prototypical ITK registration metrics, greatly increasing efficiency; (4) automated (usually hidden) parameter method. This design is overall similar to that of ITK . A few key components estimators are available; (5) transforms may include high-dimensional differ: (1) optimizers require that transforms update themselves; (2) metrics and deformation fields. One additional difference (not shown) is that “fixed” images optimizers are multi-threaded; (3) memory is shared across both optimizers and may also have a transformation, although this is not modified by the optimizer. 2.1.5. Affine and deformable similarity metrics should look as (itkQuasiNewtonOptimizerv4), multiple objective opti- similar as possible mization (itkMultiGradientOptimizerv4), or global optimization (itkMultiStartOptimizerv4). The Metricsv4 frameworksupportsthisgoalinthatitisas trivial to implement a mutual information Demons algorithm as 2.1.8. GPU and multi-core acceleration will open up new it is to implement a sum of squared differences BSpline or affine applications for image registration registration algorithm. Thus, full plug-and-play support exists See GPUPDEDeformable for a GPU example. Furthermore, across transforms. the new metric framework N cores to accelerate metric, gradi- 2.1.6. Users should be able to combine multiple similarity metrics, ent and optimization steps. A recent real-world application of the some of which may operate on different data types new Insight ToolKit implementation of the symmetric normaliza- This is achievable with the existing itkMultiGradientOpt tion algorithm showed a speed-up of almost a factor of six when imizerv4 through the multivariate itkObjectToObjec comparing single core to eight core execution time. This speed-up tMultiMetricv4 or through the multi-channel traits (itkV is achieved by multi-threading the similarity metric, the gradi- ent descent update, the regularization terms and the composition ectorImageToImageMetricTraitsv4)that allow met- rics to deal with multi-channel pixels, all of which were con- components of the method. Thus, every essential step exploits intrinsic parallelism in the algorithm. Decreased execution time tributed for ITK .The itkObjectToObjectMultiMetri cv4 was used in our winning entry of the SATA 2013 “dog leg” means more rapid turnaround for users, faster turn-around in testing and higher throughput on large-scale computing tasks. challenge. 2.1.7. Optimizers and transformations should interact flexibly 2.1.9. Improve memory efficiency in optimization framework Optimizersv4 includes optimizers that are applicable to Memory optimizations are critical for efficient use of large local both linear and deformable transformations, which include transforms. In ITK , transform parameters are no longer copied convergence monitoring and enable 2nd order optimization within the optimizer, but rather left in-place in transform. Metric Frontiers in Neuroinformatics www.frontiersin.org April 2014 | Volume 8 | Article 44 | 3 Avants et al. ITKv4 image registration FIGURE 2 | Clockwise: Define x in  and z in  as the same material Second, we remove the shape (diffeomorphic) differences (D). Consequently, I J point but existing in different domains. The point y is in a domain that is we have a composite mapping, computed via the mutual information intermediate between  and  . The standard approach in the ITKv4 similarity metric, that identifies I(x) ≈ J(A(φ(x))) = J (y) = J(z). The I J mi Affine registration framework is to map image J (B) to image I (A) by first image J (y) represents J after application of the affine transformation A Affine identifying the linear transformation, →, between the images, shown in (C). i.e., J(A(x)). Code and data for this example are here. gradient memory is shared between optimizer and metric, and begin nomenclature definition modifications by the optimizer are done in place when possible. 4 A physical point: x ∈  where  is the domain, usually of Finally, we summarize ITK changes through quantitative an image. metrics: d n An image: I :  → R where n is the number of components per pixel and d is • Over 12 new multi-threaded image registration metrics are dimensionality. A second image is J. availableinv4. Domain map: φ :  →  where → may be replaced I J with any mapping symbol below. • Four application-level registration methods, with plug-and- Affine mapping: ↔ a low-dimensional invertible play architecture, are available for high-level inclusion in transformation: affine, rigid, translation, projects such as Slicer and SimpleITK. etc. • All contributions are unit-tested and have greater than 85% Affine mapping: → designates the direction an affine code coverage, in accordance with ITK standards. mapping is applied. • A complete refactoring of the ITK transform hierarchy that Deformation field:  deformation field mapping J to I.May makes transforms thread-safe, applicable to high-dimensional not be invertible. optimization and easily used in multi-core computing. Fourty- Spline-based mapping: e.g., B-Spline field mapping J to I. one classes, in total, were impacted by this refactoring. Diffeomorphism: Represented as, these are • We added transparent vector support to two key interpola- differentiable maps with differentiable tors that are used pervasively in ITK: the nearest neighbor and inverse. Ideally, the algorithm should output the inverse and forward mapping. linear interpolators. We added two new Gaussian interpolators. Composite mapping: φ = φ (φ (x)) is defined by→ where • An example of vector support for image metrics is in 1 2 φ is of type. itkMeanSquaresImageToImageMetricv4Vecto Not invertible:  indicates a mapping that is not rRegistrationTest.cxx. invertible. Perform image warping: As an example, → J represents the Below we will discuss: (0) an organizing nomenclature matched application of an affine transform → to to the ITK framework, (1) gradient-based optimization image J such that → J = J(A(x)). within the framework, (2) techniques to estimate optimiza- Similarity measure: or ≈ indicates the metric s that s s compares a pair of images. tion parameters for arbitrary metric and transformation end nomenclature definition combinations, (3) a ITK instance implementing general- ized diffeomorphic matching, (4) several applications of the representing algorithms and applications of registration. Ideally, updated framework within different neuroinformatics-relevant any standard algorithm can be written in the nomenclature below. domains. We would then write a standard Demons registration appli- cation that maps one image, J, into the space of I (presumably 2.2. NOMENCLATURE a template) as: The nomenclature below designates an image registration algorithm symbolically. This nomenclature is intended to be a descriptive and technically consistent system for visually I → J which symbolizes I ≈ J (A (φ (x))) , Frontiers in Neuroinformatics www.frontiersin.org April 2014 | Volume 8 | Article 44 | 4 Avants et al. ITKv4 image registration with A an affine mapping and φ a generic deformation. The nota- correlation and mutual information (Hermosillo et al., 2002). tion means that the algorithm first optimizes an affine mapping, Both are implemented in ITK for transformations with local ∂J(φ(x,p)) →,between J and I. This is followed by a deformation in the sec- and global support. The term is the gradient of J at ∂φ ond stage,,from → J to I. In terms of transformation compo- ∂φ φ(x)and is the Jacobian of the transformation taken with ∂p sition, we would write → J = J (x) = J(φ (φ (x))) w Affine Demons respect to its parameter. The transform φ(x, p) may be an affine where J is the result of warping J to I.The φ are the specific map i.e., φ(x, p) = Ax + t where A is a matrix and t a transla- functions corresponding to the schematic arrows. Note, also, that tion. Alternatively, it may be a displacement field where φ(x, p) = the tail of the arrow indicates the transform’s domain. The arrow- x + u(x)and u is a vector field. In ITK , both types of maps head indicates its range. Finally, we denote the similarity metric are interchangeable and may be used in a composite transform to as ≈ which indicates a sum of squared differences (the default compute registrations that map to a template via a schematic such similarity metric). ITK supports metrics such as mutual infor- ≈ ≈ as I ≈→ J, I → J, I  → J or, mixing similarity metrics, ≈ ≈ mi b cc mation, , or cross-correlation, . We will use this nomenclature mi cc I ≈ ≈ → J . cc mi i to write schematics for registration applications in the following The most commonly used optimization algorithm for image sections. registration is gradient descent, or some variant. In the above framework, the gradient descent takes on the form of 2.3. OPTIMIZATION FRAMEWORK The general ITK optimization criterion is summarized as: ∂M ∂M φ p , x = φ p + λ , ··· , , x , new old ∂p ∂p 1 n Find mapping φ(x, p) ∈ T such that M I, J,φ(x, p) is minimized. (1) where λ is the overall learning rate and the brackets hold the vector of parameter updates. While, for functional mappings, this formulation is not strictly In addition to basic gradient descent, we implement non- linear gradient descent optimization strategies which combine the correct, the practical implementation of even high-dimensional continuous transformations involves parameterization. The space conjugate gradient or gradient descent method with line search. In ITK , we implement the classic golden section approach to T restricts the possible transformations over which to optimize the mapping φ. The arguments to φ are its parameters, p,and identifying the optimal gradient step-size at each iteration. The generic conjugate gradient approach is performed via: the spatial position, x.Notethat, in ITK , the image I may also contain a mapping, although it is not directly optimized in most 2 2 cases. As will be seen later in the document, this mapping may ∇M −∇M t t − 1 γ = , (3) also be used within large deformation metrics. ∇M t − 1 The similarity metric, M, is perhaps the most critical com- CG =∇M + γ CG , t t t − 1 ponent in image registration. Denote a parameter set as p = (p , p ... p ). The metric (or comparison function between 1 2 n where CG is the conjugate gradient. The golden section line search images) is then defined by M(I, J,φ(x, p)). For instance, M = determines the specific weight,  , of the update to the current opt I(x) − J(φ(x, p)) i.e., the sum of squared differences (SSD) parameters such that metric. Its gradient with respect to parameter p is (using the chain rule), p = p + CG  . new old t opt ∂M ∂M ∂J φ(x, p) ∂φ Note that a naive application of gradient descent will not produce M = = | . (2) p x ∂p ∂J ∂φ ∂p i i a smooth change of parameters for transformations with mixed parameter types. For instance, a change, , to parameter p will This equation provides the metric gradient specified for sum of produce a different magnitude of impact on φ if p is a translation squared differences (at point x) but similar forms arise for the rather than a rotation. Thus, we develop an estimation framework that sets “parameter scales” (in ITK parlance) which, essentially, In the Demons example above, the reader may ask: why does the affine customize the learning rate for each parameter. The update to φ mapping, A, not appear “inside” the deformable mapping, φ? Indeed, this via its gradient may also include other steps (such as Gaussian ordering of transformations is feasible. However, this is not what we typically smoothing) that project the updated transform back to space T . use in our own practice of image registration and is not what we recommend. Multi-threading is achieved in the gradient computation, trans- The reason is that we usually seek a deformable mapping into a template formation update step and (if used) the regularization by dividing space that is common across many subjects (i.e., the “tail” is in the same the parameter set into computational units that correspond to domain across subjects). This enables methods such as Jacobian-based mor- phometry and other groupwise comparison conveniences. For example, see contiguous sub-regions of the image domain. dφ Figures 1 through 4 in (Kim et al., 2008) which explains the classical approach In termsofcode, theJacobian, | , is calculated at a dp of Jacobian-based morphometry as applied to traumatic brain injury. See physical point using the function ComputeJacobianWithRespect also (Gaser et al., 2001; Riddle et al., 2004; Dubb et al., 2005; Lepore et al., ToParameters(mappedFixedPoint, Jacobian).Notethatitiseval- 2006; Rohlfing et al., 2006; Studholme and Cardenas, 2007). Additional uses uated at point x notatpoint φ(x, p). We then use the func- of Jacobian-based morphometry are shown in our examples section in the “C” example and the “asymmetry” example. tion ComputeMovingImageGradientAtPoint(mappedMovingPoint, Frontiers in Neuroinformatics www.frontiersin.org April 2014 | Volume 8 | Article 44 | 5 Avants et al. ITKv4 image registration mappedMovingImageGradient) to compute the moving image which is regularized by the linear operator L. ITK uses Gaussian gradient when there is no pre-warping. ComputeMovingImageGra smoothing which is the Green’s kernel for generalized Tikhonov dientAtPoint uses central differences (or a gradient filter) in the regularization (Nielsen et al., 1997). This objective is readily opti- dJ(φ(x,p)) mized using an approach that is similar to that proposed by Beg. moving image space to compute the image gradient, . dφ If one is doing pre-warping, then we have an index access to Generalization of the LDDMM gradient for other metrics basi- cally follows (Hermosillo et al., 2002) with a few adjustments the warped moving image. We compute the warped image J as J (x) = J(φ(x, p)). Then, to accomodate diffeomorphic mapping. Figure 3 shows an ITK result on a standard example for large deformation registra- dJ dJ φ x, p d φ(x, p) tion. We will evaluate this diffeomorphic mapping, along with = (4) parameter estimation, in a later section. dx dφ dx −1 dJ φ(x, p) dJ d φ(x, p) 2.5. PARAMETER SCALE ESTIMATION = . We choose to estimate parameter scales by analyzing the result dφ dx dx of a small parameter update on the change in the magnitude of In code, we use ComputeMovingImageGradientAtIndex(index, physical space deformation induced by the transformation. The dJ mappedMovingImageGradient) to get and transform this impact from a unit change of parameter p may be defined in dx image gradient via the inverse Jacobian by calling mappedMovingI multiple ways, such as the maximum shift of voxels or the aver- mageGradient= TransformCovariantVector(mappedMovingImag age norm of transform Jacobians (Jenkinson and Smith, 2001). eGradient, mappedMovingPoint). Denote the unscaled gradient descent update to p as p.The goal is to rescale p to q = s · p,where s is a diagonal matrix 2.4. DIFFEOMORPHIC MAPPING WITH ARBITRARY METRICS diag(s , s ... s ), such that a unit change of q will have the same 1 2 n i The framework proposed above, in general form, encom- impact on deformation for each parameter i = 1 ... n. passes both classic affine mapping as well as more recent large As an example, we want φ(x, p ) − φ(x, p )= constant new old deformation strategies. Beg proposed the Large Deformation regardless of which of the i parameters is updated by the Diffeomorphic Metric Mapping (LDDMM) algorithm (Miller unit change. The unit is an epsilon value e.g., 1.e-3. Rewrite et al., 2005) which minimizes the sum of squared differences ∂J(φ(x,p)) ∂φ ∂φ ∂M ∂M ∂M , ··· , as , ··· , . To determine the ∂p ∂p ∂J ∂φ ∂p ∂p 1 n 1 n criterion between two images. LDDMM parameterizes a diffeo- relative scale effects of each parameter, p ,wecan factor out morphism through a time varying velocity field that is integrated 4 the constant terms on the outside of the bracket. Then the through an ODE. In ITK , we implement an alternative to ∂φ modified gradient descent step becomes diag(s) . We iden- ∂p LDDMM that also uses a time varying field and an ODE but tify the values of diag(s) by explicitly computing the values of minimizes a more general objective function: φ(x, p ) − φ(x, p ) with respect to an  change. A criti- new old cal variable, practically, is which x to choose for evaluation of φ(x, p ) − φ(x, p ). The corners of the image domain work new old E(v) = M I, J,φ + w Lv  dt. (5) 1,0 t well for affine transformations. In contrast, local regions of small radius (approximately 5) work well for transformations with local This is an instance of Equation (1) where w is a scalar weight and support. Additional work is needed to verify optimal parameters φ is a standard integration of the time-varying velocity field, v , for this new ITK feature. However, a preliminary evaluation is 1,0 t FIGURE 3 | An ITK diffeomorphic mapping of the type I  J. The “C” 15 5 integration steps were used in the Runge-Kutta ODE integration over and 1/2 “C” example illustrate the large deformations that may be achieved the velocity field. A two core MacBook Air computed this registration in 110 s. with time varying velocity fields. In this case, the moving (deforming) image The images each were of size 150 × 150. See C for a reproducible example of is the 1/2 “C.” The right panels illustrate the deformed grid for the this registration and the data. In addition, we provide an example of how the transformation of the “C” to 1/2 “C” (middle right) and its inverse mapping Jacobian determinant is computed from the deformation field resulting from (far right) which takes the 1/2 “C” to the reference space. The unit time this registration via an ANTs program CreateJacobianDeterminantIm interval is discretized into 15 segments in order to compute this mapping. age. Frontiers in Neuroinformatics www.frontiersin.org April 2014 | Volume 8 | Article 44 | 6 Avants et al. ITKv4 image registration performed in the results section. The new parameter scale esti- The ANTs contribution is valuable, in part, because of mation effectively reduces the number of parameters that the user the tremendous range of registration problems that exist in must tune from k + 1(λ plus the scales for each parameter type neuroinformatics and biomedical imaging in general. While it where there are k types) to only 1, the learning rate. is not possible to solve all registration problems with a gen- The learning rate, itself, may not be intuitive for a user to eral framework, one cannot afford to invent new solutions for set. The difficulty—across problem sets—is that a good learn- every instance one encounters. Our general optimization-driven strategies have proven to be invaluable to setting performance ing rate for one problem may result in a different amount of change per iteration in another problem. Furthermore, the dis- standards in a variety of application domains. In this section, crete image gradient may become invalid beyond one voxel. Thus, we highlight some of the lesser known capabilities of ANTs and it is good practice to limit a deformation step to one voxel spac- ITK with reproducible examples that include data and specific ing (Jenkinson and Smith, 2001). We therefore provide the users commands to ANTs and/or ITK. A list of these examples follows: the ability to specify the learning rate in terms of the maxi- mum physical space change per iteration. As with the parameter 1. The Basic Brain Mapping example shows how one can map scale estimation, the domain over which this maximum change two “whole head” T1-weighted MRIs where one is a tem- is estimated impacts the outcome and similar practices are rec- plate that contains a researcher’s prior knowledge defining ommended for both cases. This feature is especially useful for the “interesting” parts of the image. Within ITK ,this allowing one to tune gradient descent parameters without being domain knowledge is used to focus the Metricsv4 on concerned about which similarity metric is being used. That is, only those parts of the image while masking the remain- it effectively rescales the term λ∂M/∂p to have a consistent effect, der. Furthermore, a second part of this example shows how for a given λ,regardlessofthe metric choice.Incombination with the ITK composite transform may be used to initialize new our non-linear conjugate gradient approach (our current opti- registration solutions as well as how masking functionality mization of choice for linear registration), this strategy drastically may be employed to ignore information that is irrelevant or reduces the parameter setting burden for users. obstructive to the registration optimization. We have previ- ously employed these strategies in brain mapping with lesions 3. RESULTS (Avants et al., 2008; Kim et al., 2008, 2010, 2013; Tustison et 3.1. EXAMPLE APPLICATIONS OF THE ITK FRAMEWORK al., 2011). As part of our work in ITK refactoring, we built, in paral- 2. We use updated ITK methods in template construction with lel to library programming, an application interface that allows a reproducible example based on face and brain data: ANTs high-level access to the deep layers of ITK registration. These cur- template construction. This work has been employed in rently exist in the Advanced Normalization Tools (ANTs) software different species, age-ranges, and imaging modalities. The (link). While ANTs still serves as intermediate (vs. direct) access- resulting template is an image that captures the expected point to these tools, it provides a high-degree of customization shape and appearance as defined by the population sample, possibilities simply through a command line interface and script- transformation model and intensity comparison metric. ing. Therefore, a user is not required to write new low-level 3. A large deformation example implementing the classic “letter (interestingly, C++ is now considered “low-level”) software. C” example provided, originally, by Gary Christensen. While Despite its relative youth, the ANTs wrapping of ITK func- extremely flexible, these algorithms have not found a unique tionality has been employed with notable success in recent public, identity in terms of translational applications yet remain of unbiased, international evaluation studies. ANTs was instrumen- theoretical interest. This example shows a user how to define tal to a first-place finish in SATA 2013 in two of three categories the parameters of a registration based on optimizing a time (based on the median performance) where the ANTs approach varying velocity field. was considerably simpler than that employed by close finishers. 4. We present a separate example of how to compute landmark- While the evaluation of deep-gray matter registration showed based registration error. ITK uses LPS coordinates to rep- relatively subtle differences, the ANTs solution to the multivari- resent physical space. If you need to convert landmarks to ate canine leg data outclassed all other entrants. Notably, the physical space, see the discussion here: LPS physical space. ANTs solution used a multiple metric approach that simultane- We have an example illustrating how to change point coor- ously compared two modalities during registration as in Avants dinates and apply ITK transforms to landmarks here. This et al. (2008). In the cardiac data, the ANTs solution was the exercise can be useful for landmark-based registration or in only one that was fully automated resulting in a ≈15% per- evaluating registration accuracy. formance loss which can easily be overcome by a modicum of 5. We show how to perform motion correction to time series user intervention. Furthermore, ANTs/ITK-based methods fin- data here although we do not claim this approach is opti- ished a clear first-place in the BRATS 2013 challenge. Our entry mal. The method registers each frame from a 4D time series used intensity asymmetry as a key feature to segment brain to a fixed reference image and stores the resampled set in a tumors based on multiple modality MRI. Thus, these methods new 4D image. All transformation parameters are stored in a are within the leading ranks of image registration methodolo- corresponding csv file. gies as evaluated in recent work as well as in the more traditional 6. An advanced example with heavy use of statistics via R and brain (Klein et al., 2010) and lung CT (Murphy et al., 2011) ANTsR is in a study of public test-retest fmri data. This study studies. is not published and may be subject to change. Frontiers in Neuroinformatics www.frontiersin.org April 2014 | Volume 8 | Article 44 | 7 Avants et al. ITKv4 image registration 7. The classic car example shown in ANTs talks is here. This 3.2.1. Parameter estimation across metrics illustrates the benefit of mutual information in deformably ITK provides similarity metrics that may be applied for both mapping wild-type images and highlights the fact that deformable and affine registration. In a previous section, we pro- ITK applications exist outside of medical imaging. vided a parameter estimation strategy that is applicable to both 8. A basic multistart example. A more advanced example deformable and affine transformations with arbitrary metrics. for brain mapping with a template mask is available in Denote images I, J, K, where the latter two are “moving” images, antsCorticalThickness.sh. These optimization methods over- and K is an intensity-inverted version of J.Wethenevaluatethe come local minima by running registration searches from following schema, a variety of starting points and greedily storing the best solution. I ≈→ J, I ≈ → K, I ≈ → K cc mi 9. This asymmetry analysis example uses the “virtual domain” where, for each schematic, we use the corresponding metric for feature to reduce bias caused by mapping an image asymmet- both affine and diffeomorphic mapping. Furthermore, we keep rically to a reference. Note that we measure the point-wise the same parameters for each registration by exploiting parameter asymmetry, in this example, via the Jacobian determinant scale estimators. Figure 5 shows the candidate images for this test. image as in Kim et al. (2008). If one repeats this analy- As shown in Figure 5, very similar results are achieved for each sis across a population—and maps the Jacobian measure- schematic without additional parameter tuning. To determine ments of asymmetry to a common space—then one may this quantitatively, we perform registration for each schematic perform a statistical analysis of population-level asymme- and then compare the Dice overlap of a ground-truth three-tissue try. Longitudinal analysis and asymmetry analysis potentially segmentation. For each result, we have the Dice overlap of dark suffer from the same confound Yushkevich et al. (2010). A tissue (cerebrospinal fluid, CSF), medium intensity tissue (gray related longitudinal mapping script is here. matter) and bright tissue (white matter). For the mean squares 10. We also show how manual labelings can be used to restrict metric, we have: 0.588, 0.816, and 0.90; for CC, we have: 0.624, registration in a challenging registration scenario (this exam- 0.786, 0.882; for MI, we have: 0.645, 0.779, 0.858. Mutual infor- ple will be improved in the future): registration guided by mation does best for the CSF while mean squares does best for (crude) labels. other tissues. CC performs in the mid-range for all classes of 11. A simple orange to apple RGB image registration example for tissue. Thus, a single set of tuned parameters provides a reason- color images is listed at: itkMeanSquaresImageToImageMet- able result for an affine plus diffeomorphic mapping across three ricv4VectorRegistrationTest. If one compiles the ITK tests, different metrics. While improvement might be gained by fur- then this example can be run to produce Figure 4. ther tuning for each metric, this result shows that our parameter estimation method achieves the goal of reducing user burden. These examples cover many applications for which no “ground 3.2.2. Automated Brain Labeling Task truth” evaluation data exists. The next section seeks to add some quantitative reference to these examples. First, we show flexibility All R and bash analysis scripts for this section are here: https://github.com/stnava/ITKv4Documentation/tree/frontiers/sc and consistency of our framework in a simple example comparing registration with a variety of metrics and a consistent parameter ripts. The ITK core functionality formed the heart of the ref- erence results provided for the SATA2013 challenge at MICCAI set. Second, we quantify the benefit of ITK registration in com- parison to a method implemented based upon ITK registration 2013. In this sense, these methods have been heavily evaluated on both basic brain mapping challenges (SATA2013’s dien- technology. cephalon challenge in which ITK -based methods finished first), multivariate registration challenges (the canine MRI / dog 3.2. EVALUATION leg challenge of SATA2013 in which ITK -based methods were We first investigate the ability of our automated parameter esti- overwhelmingly the top finisher) and in the cardiac challenge mation to facilitate parameter tuning across metrics. Second, we 4 3 (in which ITK -based methods were the only fully automated compare ITK and ITK registration implementations with approach). However, for completeness, we provide an additional respect to a standard automated brain labeling task. evaluation here which focuses on comparison to a ITK method, BRAINSFit, in a different dataset than previously used to evaluate ANTs or ITK . As ground truth, we use T1 MRI data from 33 2-year old subjects as described in Gousias et al. (2008) and available at http://www.brain-development.org. Each subject’s brain is man- ually parcellated into 83 distinct regions that include ventricles, cortical areas, white matter and deep gray matter regions such as the amygdala, hippocampus and thalamus. One benefit of this data is that some of these anatomical regions are rela- FIGURE 4 | This RGB image registration example employs ITK code tively easy to align (the caudate) whereas others are relatively that repurposes a scalar metric (itkMeanSquaresImageToImageMetr difficult to align due to their small size (amygdala) or incon- icv4) for multichannel registration. sistent shape across subjects (the inferior frontal gyrus). Thus, Frontiers in Neuroinformatics www.frontiersin.org April 2014 | Volume 8 | Article 44 | 8 Avants et al. ITKv4 image registration FIGURE 5 | Three reference images, I (left), J (middle top), and K (right affine and deformable case. Thus, two parameters had to be optimized. We top), are used to illustrate the robustness of our parameter scale then applied these same parameters to register I and K via both correlation estimation for setting consistent parameters across both metrics and and mutual information. The resulting registrations (bottom row) were all of transform types. K is the negation of J andisusedtotestthe correlation similar quality. Further, the same metric is used for both affine and and mutual information registrations. We optimized, by hand, the step-length diffeomorphic mapping by exploiting the general optimization process given parameters for one metric (the sum of squared differences) for both the in Equation (1). FIGURE 6 | We compare a ITK composite schema as I ≈ ≈ → J example ANTs-based large-deformation result from the dataset is shown for cc mi i for mapping a set of {J } images to a template I to a ITK schema: illustration where we render the extracted brains as well as show select axial I ≈  ≈ → J . We use this schematic in a registration-based slices. All registrations were run on the original MRI data with no mi b mi i segmentation of multiple brain structures in a pediatric population as a preprocessing except what is done by ANTs or BRAINSFit internally. Overlap benchmark for algorithm performance, similar to Klein et al. (2010). An improvement from v3 to v4, quantified via paired t-test, is highly significant. Frontiers in Neuroinformatics www.frontiersin.org April 2014 | Volume 8 | Article 44 | 9 Avants et al. ITKv4 image registration we anticipate that performance gains due to new technology in contained within the ANTs repository. The git hashtag for the ITK will be most prominent in the more variable and chal- BRAINSFit version https://github.com/BRAINSia/BRAINSTools lenging regions. Figure 6 summarizes the study nomenclature used in this evaluation is ad7e114ab1c92bd800819b80e054825 and shows a single image pair selected from this data along 9398931c8. Both programs were run on a MacBook Pro running with the registration result given by ITK . Figure 7 summa- OS X 10.9 (13A3028) with a 2.6 GHz Intel Core i7 and 16 GB of rizes these evaluation results. The scripts for running this study RAM. are available at https://github.com/stnava/ITKv4Documentation/ To help isolate which subject pairs to deformably register, tree/frontiers/scripts. The git hashtag for the ANTs version used we first clustered the initial dataset based on an all pairs affine in this evaluation https://github.com/stnava/ANTs is ce8b5a741 registration which revealed five representative subject clusters. 4ae9e389071d756c5f36ee6cecbcfd8. The associated ITK tag is The subject pairings used during evaluation were chosen such FIGURE 7 | Above, a barplot shows the mean Dice score for each region ITK SyN algorithm, with its classic neighborhood correlation metric, and each algorithm, sorted by ANTs performance. Below, we use star outperforms BRAINSFit in several regions and more strongly in some subject plots of per-brain-region Dice overlap to compare, for each subject, the pairs than others. The legend for the plots is at lower right and shows the 4 3 ITK implementation of SyN with the ITK -based BRAINSFit algorithm. The maximum possible value for each region. Frontiers in Neuroinformatics www.frontiersin.org April 2014 | Volume 8 | Article 44 | 10 Avants et al. ITKv4 image registration Table 1 | Dice overlap for ANTs and BRAINSFit where only regions that each subject pair contained the most representative subject with q-value < 0.01 are shown. for one cluster paired with the most representative subject from another cluster; thus, the study design allows us to focus, in a meanants meanbfit FDR-adjusted principled manner, on a set of representative shape comparisons p-value where the comparisons are made across different image types. The new methods in ITK show enhanced performance within R Hippocampus 0.835 0.770 0.0083 all registration pairs. The mean overall performance gain was L Hippocampus 0.834 0.782 0.0063 approximately 6.3% with a standard deviation of 5% and T- R Amygdala 0.843 0.752 0.0027 statistic/p-value, over all structures, of 12.6 / p < 1.e − 16. We L Amygdala 0.851 0.773 0.0037 also identified which regions were most improved in ITK vs. L AnteriorTemporalLobeMe 0.849 0.752 0.0005 ITK . These regions include the left and right insula, the brain- dialPart stem, the superior temporal gyrus, parahippocampal gyrus, puta- R AnteriorTemporalLobeLate 0.754 0.633 0.0002 men, and the substantia nigra. Table 1 lists all structures and ralPart the mean Dice score for each algorithm, along with the p-value. L AnteriorTemporalLobeLat 0.786 0.668 0.0083 Figure 7 summarizes all of these findings by using star plots to eralPart R Gyri parahippocampalis et 0.813 0.756 0.0037 visualize the Dice results for every region in every subject. ambiens We also recorded the amount of processing time spent on each L Gyri parahippocampalis et 0.823 0.749 0.0007 subject, for each algorithm. Noting that the ITK algorithm also ambiens provides a dense and high-resolution forward and inverse trans- R Superior temporal gyrus 0.892 0.827 0.0001 form and does explicit transformation regularization to guarantee posterior a diffeomorphism, the algorithm takes, on average, 5 times as long 3 R Medial and inferior 0.891 0.823 0.0007 as the ITK BRAINSFit algorithm (≈10 min), assuming default temporal gyri settings. Much of this time, in ANTs, is taken up by full resolution R Lateral occipitotemporal 0.801 0.690 0.0023 image registration. If this fine level is avoided, then the disparity gyrus (gyrus fusiformis) in timing reduces to less than a factor of two, without much loss R Cerebellum 0.953 0.937 0.0066 in accuracy. Note also that ANTs and BRAINSFit each use a differ- Brainstem 0.944 0.923 0.0001 ent multithreading strategies, similarity metric implementations, R Insula 0.910 0.864 0.0001 rigid/affine registration mechanisms and optimizers making this L Insula 0.915 0.863 0.0001 overall comparison less than ideal. L Occipital lobe 0.895 0.819 0.0049 L Posterior temporal lobe 0.912 0.862 0.0001 4. DISCUSSION L Parietal lobe 0.885 0.815 0.0052 ITK is a community built and maintained toolkit and is a public R Putamen 0.909 0.869 0.0002 resource for reproducible methods. The updated ITK registra- L Putamen 0.911 0.862 0.0003 tion framework provides a novel set of user-friendly parameter R Pallidum 0.838 0.721 0.0037 setting tools and benchmark implementations of both standard L Pallidum 0.846 0.766 0.0049 and advanced algorithms. Robustness with respect to parameter L Lingual gyrus 0.876 0.795 0.0037 settings has long been a goal of image registration and ITK takes L Cuneus 0.825 0.684 0.0083 valuable steps toward the direction of automated parameter L Lateral orbital gyrus 0.728 0.598 0.0037 selection. The primary decision left up to the user, currently, L Substantial nigra 0.831 0.741 0.0023 is the feature scale at which registration should be performed. L Superior temporal gyrus 0.818 0.752 0.0043 E.g., whether the registration should focus on coarse features, anterior part fine features, etc and the different resolutions at which this should be done. While we have provided a reproducible reference comparison of registration-based brain labeling in this paper, we implementations for speed and memory. In time, it may be intend to have a more extensive series of benchmark performance possible to extend the design philosophy used here to GPU studies completed on datasets beyond the brain. However, the implementations. However, our ability to interface low and number of possible applications exceeds what can possibly be high-dimensional transformations depends heavily on generic evaluated by the core ITK developer community. Community programming. This style is less well-developed (and less well involvement is needed in order to increase the number of possible understood) in GPU applications which depend, to some registration applications and metric/transform/optimizer/data extent, on specialization. The current framework is amenable combinations that have been evaluated. At the same time, to groupwise registration strategies when used in combina- documentation, usability and examples must be pro- tion with a computing cluster. However, single core group- vided by the development team in order to improve user wise strategies are not currently implemented although one involvement. may consider basing an implementation on exisiting multi- 4.1. FUTURE WORK metric/multivariate registration tools within the current code Future work will enhance the depth and breadth of docu- base. While ITK does contain a statistics infrastructure, we mentation as well as seek to further optimize the current currently prefer using R and ANTsR for analyzing our data. Frontiers in Neuroinformatics www.frontiersin.org April 2014 | Volume 8 | Article 44 | 11 Avants et al. ITKv4 image registration However, the lack of visualization methods in ITK means volumetry of brain lateral ventricles in mri. Neuroimage 13(6 pt 1), 1140–1145. doi: 10.1006/nimg.2001.0771 that one must still move to another package to look at Gousias, I. S., Rueckert, D., Heckemann, R. A., Dyet, L. E., Boardman, one’s results. Therefore, direct interfaces to R remain useful. J. P., Edwards, A. D., et al. (2008). Automatic segmentation of brain mris SimpleITK also has a promising R interface that is similar to of 2-year-olds into 83 regions of interest. Neuroimage 40, 672–684. doi: ANTsR. 10.1016/j.neuroimage.2007.11.034 A primary challenge to the future of ITK includes, beyond Hermosillo, G., Chefd’Hotel, C., and Faugeras, O. (2002). A variational approach to multi-modal image matching. Intl.J.Comp. Vis. 50, 329–343. doi: documentation, reduced C ++ fluency. As ITK leverages 10.1023/A:1020830525823 several advanced features of C ++, even experienced devel- Jack, C. R. Jr., Knopman, D. S., Jagust, W. J., Shaw, L. M., Aisen, P. S., opers may find it difficult to contribute meaningfully to the Weiner, M. W., et al. (2010). Hypothetical model of dynamic biomark- ITK software base. Therefore, the ITK community must also ers of the Alzheimer’s pathological cascade. Lancet Neurol. 9, 119–128. doi: 10.1016/S1474-4422(09)70299-6 seek to educate potential future contributors not only on ITK Jenkinson, M., and Smith, S. (2001). A global optimisation method for robust affine but also, at times, on the fundamentals or advanced exten- registration of brain images. Med. Image Anal. 5, 143–156. doi: 10.1016/S1361- sions of C ++. A second major hurdle is that ITK includes 8415(01)00036-6 a host of generic registration ingredients. However, many of Kikinis, R., and Pieper, S. (2011). 3d slicer as a tool for interactive brain tumor the most compelling new application domains require special- segmentation. Conf. Proc. IEEE Eng. Med. Biol. Soc. 2011, 6982–6984. doi: 10.1109/IEMBS.2011.6091765 ization. Specialization may be needed for a specific imaging Kim, J., Avants, B., Patel, S., Whyte, J., Coslett, B. H., Pluta, J., et al. (2008). modality, via hardware interface or in the use of domain- Structural consequences of diffuse traumatic brain injury: a large defor- specific prior knowledge. Therefore, we envision the next phase mation tensor-based morphometry study. Neuroimage 39, 1014-1026. doi: of ITK development may focus on using the toolkit to sup- 10.1016/j.neuroimage.2007.10.005 port its specialization in solving high-impact and translational Kim, J., Avants, B., Whyte, J., and Gee, J. C. (2013). Methodological considerations in longitudinal morphometry of traumatic brain injury. Front. Hum. Neurosci. applications. Hopefully, this transition will occur in the near 7:52. doi: 10.3389/fnhum.2013.00052 future. Kim, J., Whyte, J., Patel, S., Avants, B., Europa, E., Wang, J., et al. (2010). Resting cerebral blood flow alterations in chronic traumatic brain injury: an arte- rial spin labeling perfusion fmri study. J. Neurotrauma 27, 1399-1411. doi: REFERENCES 10.1089/neu.2009.1215 Ackerman, M. J., and Yoo, T. S. (2003). The visible human data sets (VHD) and insight toolkit (ITk): experiments in open source software. AMIA Annu. Symp. Klein, S., Staring, M., Murphy, K., Viergever, M. A., and Pluim, J. P. W. (2010). elastix: a toolbox for intensity-based medical image registra- Proc. 2003:773. Avants, B. B., Tustison, N. J., Song, G., Cook, P. A., Klein, A., and Gee, J. C. (2011). tion. IEEE Trans. Med. Imaging 29, 196–205. doi: 10.1109/TMI.2009. A reproducible evaluation of ANTs similarity metric performance in brain image registration. Neuroimage 54, 2033–2044. doi: 10.1016/j.neuroimage. Lepore, N., Brun, C. A., Chiang, M.-C., Chou, Y.-Y., Dutton, R. A., Hayashi, K. M., et al. (2006). Multivariate statistics of the jacobian matrices in tensor based mor- 2010.09.025 phometry and their application to hiv/aids. Med. Image Comput. Comput. Assist. Avants, B., Duda, J. T., Kim, J., Zhang, H., Pluta, J., Gee, J. C., et al. (2008). Multivariate analysis of structural and diffusion imaging in traumatic brain Interv. 9(pt 1), 191–198. Metz, C. T., Klein, S., Schaap, M., van Walsum, T., and Niessen, W. J. (2011). injury. Acad. Radiol. 15, 1360-1375. doi: 10.1016/j.acra.2008.07.007 Baloch, S., and Davatzikos, C. (2009). Morphological appearance mani- Nonrigid registration of dynamic medical imaging data using nd + t b-splines and a groupwise optimization approach. Med. Image Anal. 15, 238–249. doi: folds in computational anatomy: groupwise registration and morphologi- cal analysis. Neuroimage 45(1 Suppl), S73–S85. doi: 10.1016/j.neuroimage. 10.1016/j.media.2010.10.003 Miller, M. I., Beg, M. F., Ceritoglu, C., and Stark, C. (2005). Increasing the power 2008.10.048 Bearden, C. E., van Erp, T. G. M., Dutton, R. A., Tran, H., Zimmermann, L., Sun, of functional maps of the medial temporal lobe by using large deformation dif- feomorphic metric mapping. Proc. Natl. Acad. Sci. U.S.A. 102, 9685–9690. doi: D., et al. (2007). Mapping cortical thickness in children with 22q11.2 deletions. Cereb. Cortex 17, 1889–1898. doi: 10.1093/cercor/bhl097 10.1073/pnas.0503892102 Murphy, K., van Ginneken, B., Reinhardt, J. M., Kabus, S., Ding, K., Deng, X., Buckheit, J. B., and Donoho, D. L. (1995). “WaveLab and reproducible research,” in Wavelets and statistics,eds A. Antoniadis andG.Oppenheim (New York,NY: et al. (2011). Evaluation of registration methods on thoracic CT: the EMPIRE10 challenge. IEEE Trans. Med. Imaging 30, 1901–1920. doi: 10.1109/TMI.2011. Springer), 55–81. Chen, M., Lu, W., Chen, Q., Ruchala, K. J., and Olivera, G. H. (2008). A simple 2158349 fixed-point approach to invert a deformation field. Med. Phys. 35, 81–88. doi: Nielsen, M., Florack, L., and Deriche, R. (1997). Regularization, scale- space, and edge detection filters. J. Math. Imaging Vis. 7, 291–307. doi: 10.1118/1.2816107 Cheung, M. R., and Krishnan, K. (2009). Interactive deformation registration of 10.1023/A:1008282127190 Peyrat, J.-M., Delingette, H., Sermesant, M., Xu, C., and Ayache, N. (2010). endorectal prostate mri using itk thin plate splines. Acad. Radiol. 16, 351–357. doi: 10.1016/j.acra.2008.09.011 Registration of 4d cardiac ct sequences under trajectory constraints with multi- channel diffeomorphic demons. IEEE Trans. Med. Imaging 29, 1351–1368. doi: Christensen, G. E., Rabbitt, R. D., and Miller, M. I. (1996). Deformable templates using large deformation kinematics. IEEE Trans. Image Process. 5, 1435–1447. 10.1109/TMI.2009.2038908 Riddle, W. R., Li, R., Fitzpatrick, J. M., DonLevy, S. C., Dawant, B. M., doi: 10.1109/83.536892 Dubb, A., Xie, Z., Gur, R., Gur, R., and Gee, J. (2005). Characterization of brain and Price, R. R. (2004). Characterizing changes in mr images with color- coded jacobians. Magn. Reson. Imaging 22, 769–777. doi: 10.1016/j.mri. plasticity in schizophrenia using template deformation. Acad. Radiol. 12, 3–9. doi: 10.1016/j.acra.2004.06.009 2004.01.078 Rohlfing, T., Sullivan, E. V., and Pfefferbaum, A. (2006). Deformation-based Fedorov, A., Li, X., Pohl, K. M., Bouix, S., Styner, M., Addicott, M., et al. (2011). Atlas-guided segmentation of vervet monkey brain MRI. Open Neuroimag. J. 5, brain morphometry to track the course of alcoholism: differences between intra-subject and inter-subject analysis. Psychiatry Res. 146, 157–170. doi: 186–197. doi: 10.2174/1874440001105010186 10.1016/j.pscychresns.2005.12.002 Floca, R., and Dickhaus, H. (2007). A flexible registration and evalua- tion engine (f.r.e.e.). Comput. Methods Programs Biomed. 87, 81–92. doi: Rueckert, D., Sonoda, L. I., Hayes, C., Hill, D. L., Leach, M. O., and Hawkes, D. J. (1999). Nonrigid registration using free-form deformations: applica- 10.1016/j.cmpb.2007.04.009 Gaser, C., Nenadic, I., Buchsbaum, B. R., Hazlett, E. A., and Buchsbaum, M. S. tion to breast mr images. IEEE Trans. Med. Imaging 18, 712–721. doi: 10.1109/42.796284 (2001). Deformation-based morphometry and its relation to conventional Frontiers in Neuroinformatics www.frontiersin.org April 2014 | Volume 8 | Article 44 | 12 Avants et al. ITKv4 image registration Shelton, D., Stetten, G., Aylward, S., Ibáñez, L., Cois, A., and Stewart, Yushkevich, P. A., Avants, B. B., Das, S. R., Pluta, J., Altinay, M., Craige, C. (2005). Teaching medical image analysis with the insight C., et al. (2010). Bias in estimation of hippocampal atrophy using toolkit. Med. Image Anal. 9, 605–611. doi: 10.1016/j.media.2005. deformation-based morphometry arises from asymmetric global normal- 04.011 ization: an illustration in adni 3 t mri data. Neuroimage 50, 434-445. doi: Studholme, C., and Cardenas, V. (2007). Population based analysis of directional 10.1016/j.neuroimage.2009.12.007 information in serial deformation tensor morphometry. Med. Image Comput. Comput. Assist. Interv. 10(pt 2), 311–318. Conflict of Interest Statement: The authors declare that the research was con- Tustison, N., Avants, B., Cook, P., Kim, J., Whyte, J., Gee, J., et al. (2011). ducted in the absence of any commercial or financial relationships that could be “Multivariate analysis of diffusion tensor imaging and cortical thickness maps construed as a potential conflict of interest. in a traumatic brain injury (tbi) cohort using advanced normalization tools (ants) Journal of Neurotrauma 28 A111,” in 29th Annual National Neurotrauma Received: 01 November 2013; accepted: 30 March 2014; published online: 28 April Symposium (Hollywood Beach, FL). 2014. van Dalen, J. A., Vogel, W., Huisman, H. J., Oyen, W. J. G., Jager, G. J., Citation: Avants BB, Tustison NJ, Stauffer M, Song G, Wu B and Gee JC (2014) The and Karssemeijer, N. (2004). Accuracy of rigid CT-FDG-PET image registra- Insight ToolKit image registration framework. Front. Neuroinform. 8:44. doi: 10.3389/ tion of the liver. Phys. Med. Biol. 49, 5393–5405. doi: 10.1088/0031-9155/49/ fninf.2014.00044 23/014 This article was submitted to the journal Frontiers in Neuroinformatics. Vercauteren, T., Pennec, X., Perchant, A., and Ayache, N. (2009). Diffeomorphic Copyright © 2014 Avants, Tustison, Stauffer, Song, Wu and Gee. This is an open- demons: efficient non-parametric image registration. Neuroimage 45(1 Suppl), access article distributed under the terms of the Creative Commons Attribution S61–S72. doi: 10.1016/j.neuroimage.2008.10.040 License (CC BY). The use, distribution or reproduction in other forums is permit- Yoo, T. S., Ackerman, M. J., Lorensen, W. E., Schroeder, W., Chalana, V., Aylward, ted, provided the original author(s) or licensor are credited and that the original S., et al. (2002). Engineering and algorithm design for an image processing api: publication in this journal is cited, in accordance with accepted academic prac- a technical report on itk–the insight toolkit. Stud. Health Technol. Inform. 85, tice. No use, distribution or reproduction is permitted which does not comply with 586–592. these terms. Frontiers in Neuroinformatics www.frontiersin.org April 2014 | Volume 8 | Article 44 | 13 http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Frontiers in Neuroinformatics Pubmed Central

The Insight ToolKit image registration framework

Frontiers in Neuroinformatics , Volume 8 – Apr 28, 2014

Loading next page...
 
/lp/pubmed-central/the-insight-toolkit-image-registration-framework-l4J085pjzO

References (87)

Publisher
Pubmed Central
Copyright
Copyright © 2014 Avants, Tustison, Stauffer, Song, Wu and Gee.
ISSN
1662-5196
eISSN
1662-5196
DOI
10.3389/fninf.2014.00044
Publisher site
See Article on Publisher Site

Abstract

ORIGINAL RESEARCH ARTICLE published: 28 April 2014 NEUROINFORMATICS doi: 10.3389/fninf.2014.00044 1 2 1 1 1 Brian B. Avants *, Nicholas J. Tustison , Michael Stauffer , Gang Song , Baohua Wu and James C. Gee Penn Image Computing and Science Laboratory, Department of Radiology, University of Pennsylvania, Philadelphia, PA, USA Department of Radiology and Medical Imaging, University of Virginia, Charlottesville, VA, USA Edited by: Publicly available scientific resources help establish evaluation standards, provide a Hans J. Johnson, The University of platform for teaching and improve reproducibility. Version 4 of the Insight ToolKit ( ITK ) Iowa, USA seeks to establish new standards in publicly available image registration methodology. Reviewed by: 4 4 ITK makes several advances in comparison to previous versions of ITK. ITK supports Marcel Prastawa, University of both multivariate images and objective functions; it also unifies high-dimensional Utah, USA Andrey Fedorov, Brigham and (deformation field) and low-dimensional (affine) transformations with metrics that are Women’s Hospital, USA reusable across transform types and with composite transforms that allow arbitrary *Correspondence: series of geometric mappings to be chained together seamlessly. Metrics and optimizers Brian B. Avants, Department of take advantage of multi-core resources, when available. Furthermore, ITK reduces the Radiology, University of parameter optimization burden via principled heuristics that automatically set scaling Pennsylvania, 3600 Market St., Suite 370, Philadelphia, PA 19104, USA across disparate parameter types (rotations vs. translations). A related approach also e-mail: avants@grasp.cis.upenn.edu constrains steps sizes for gradient-based optimizers. The result is that tuning for different metrics and/or image pairs is rarely necessary allowing the researcher to more easily focus on design/comparison of registration strategies. In total, the ITK contribution is intended as a structure to support reproducible research practices, will provide a more extensive foundation against which to evaluate new work in image registration and also enable application level programmers a broad suite of tools on which to build. Finally, we contextualize this work with a reference registration evaluation study with application to pediatric brain labeling. Keywords: registration, MRI, brain, open-source, death 1. INTRODUCTION resources have become critical to setting performance standards in international challenges that evaluate “real world” registra- As image registration methods mature—and their capabilities tion scenarios (see, for instance, the SATA 2013 and BRATS 2013 become more widely recognized—the number of applications challenges at MICCAI in Nagoya, Japan). increase (Rueckert et al., 1999; van Dalen et al., 2004; Miller et al., 2005; Shelton et al., 2005; Chen et al., 2008; Baloch and One source of benchmark methodology is the Insight ToolKit (ITK) (Yoo et al., 2002; Ackerman and Yoo, 2003), which marked Davatzikos, 2009; Cheung and Krishnan, 2009; Peyrat et al., 2010; Fedorov et al., 2011; Kikinis and Pieper, 2011; Metz a significant contribution to medical image processing when it et al., 2011; Murphy et al., 2011). Consequently, image regis- first emerged at the turn of the millennium. Since that time, ITK tration transitioned from being a field of active research, and has become a standard-bearer for image processing algorithms and, in particular, for image registration methods. In a review few applied results, to a field where the main focus is trans- lational. Image registration is now used to derive quantita- of ITK user interests, image registration was cited as the most important contribution of ITK (personal communication with tive biomarkers from images (Jack et al., 2010), plays a major role in business models and clinical products (especially in Terry Yoo). Numerous papers use ITK algorithms as standard ref- erences for implementations of Demons registration and mutual radiation oncology) (Cheung and Krishnan, 2009), has led to numerous new findings in studies of brain and behavior (e.g., information-based affine or B-Spline registration (van Dalen Bearden et al., 2007) and is a critical component in applica- et al., 2004; Shelton et al., 2005; Floca and Dickhaus, 2007; Chen tions in pathology, microscopy, surgical planning, and more et al., 2008; Cheung and Krishnan, 2009). Multiple toolkits extend (Miller et al., 2005; Shelton et al., 2005; Floca and Dickhaus, ITK registration methods in unique ways. Elastix provides very 2007; Chen et al., 2008; Cheung and Krishnan, 2009; Peyrat fast and accurate B-Spline registration (Klein et al., 2010; Murphy et al., 2010; Kikinis and Pieper, 2011; Murphy et al., 2011). et al., 2011). The diffeomorphic demons is a fast/efficient approx- Despite the increasing relevance of image registration across imation to a diffeomorphic mapping (Vercauteren et al., 2009). application domains, there are relatively few reference algorithm ANTs provides both flexibility and high average performance (Avants et al., 2011). The BRAINSFit algorithm is integrated into implementations available to the community. Furthermore, these Slicer for user-guided registration (Kikinis and Pieper, 2011). Each of these toolkits has both strengths and weaknesses (Klein This work is supported by National Library of Medicine sponsored ARRA stimulus funding. et al., 2010; Murphy et al., 2011) and was enabled by an ITK core. Frontiers in Neuroinformatics www.frontiersin.org April 2014 | Volume 8 | Article 44 | 1 Avants et al. ITKv4 image registration The Insight ToolKit began a major refactoring effort in 2010. optimization technique (in the Optimizersv4 directory), the The refactoring aimed to both simplify and extend the tech- metric measuring the registration quality (the Metricsv4 niques available in version 3.x with methods and ideas from a directory), the images or other data objects that enter the met- new set of prior work (Christensen et al., 1996; Rueckert et al., ric and the parameters that are being optimized. The parame- 1999; Jenkinson and Smith, 2001; Miller et al., 2005; Peyrat et al., ters are usually defined by a geometric transformation but may 2010; Avants et al., 2011). To make this technology more accessi- point to other relevant objects. Any of ITK ’s transforma- ble, ITK unifies the dense registration framework (displacement tions may be optimized by the framework. New transformations, field, diffeomorphisms) with the low-dimensional (B-Spline, relative to ITK ,include the DisplacementField trans- affine, rigid) framework by introducing composite transforms, forms that are useful for engendering Demons or B-Spline reg- deformation field transforms, and specializations that allowed istration strategies. New VelocityField transforms are also these to be optimized efficiently. A sub-goal set for ITK was available. A typical application developer would employ all of to simplify parameter setting by adding helper methods that these components. A good starting point for new users who use well-known principles of image registration to automatically wish to see how these tools work together, in source code, is scale transform components and set optimization parameters. foundinthe tests. Forinstance,see thefiles itkTimeVaryin ITK transforms are also newly applicable to objects such as vec- gBSplineVelocityFieldImageRegistrationTest. tors and tensors and will take into account covariant geometry cxx for an example of a B-Spline diffeomorphism applica- if necessary. Finally, ITK reconfigures the registration frame- tion, itkSyNImageRegistrationTest.cxx to see SyN in work to maximize multi-threading resources where possible. The ITK and itkSimpleImageRegistrationTest2.cxx revised registration framework within ITK is more thoroughly for a more basic example. integrated across transform models, is thread-safe and provides Several usability goals spurred ITK development. We sum- broader functionalitythaninprior releases. marize these here. David Donoho once commented (in paraphrase) that aca- demic publications are merely “advertisements” for the real work 2.1.2. Image registration should be achievable in one step which is constituted by the “complete instruction set” that pro- This overarching goal is best illustrated by Registration duces the results reported in the publication (Buckheit and Methodsv4 in which a user may string together a series Donoho, 1995). The first part of the remainder of this docu- of registration tools to perform (for instance) a transla- ment will provide an “advertisement” for the ITK framework and tion registration, followed by an affine registration, followed 3 4 summarize its evolution from ITK to ITK . We then detail by a diffeomorphic mapping each of which might use a potential applications of this ITK framework in the context different image similarity metric. The different transforms of a general nomenclature. While this work is indeed incom- are accumulated in the new itkCompositeTransform plete, in the sense of Donoho, we refer to source code and data which chains transforms together as in Figure 2. Thus, this when relevant. Furthermore, section 3.1 shows a series of repro- framework provides unprecedented ability to perform com- ducible examples of ITK in action. Several areas relevant to plex and staged registration mappings. Furthermore, the neuroinformatics are highlighted in these examples: optimal tem- frameworks automated parameter scaling, itkRegistrat plate construction, “challenging” registration scenarios involving ionParameterScalesEstimator, vastly reduces the dif- brain mapping in the presence of lesions or resection, registra- ficulty of tuning parameters for different transform/metric tion when initialization priors are weak, asymmetry analyses, combinations. functional MRI, and non-traditional registration strategies are all highlighted. We also establish performance benchmarks for the 2.1.3. ITK Transforms should be unified current ITK registration, in comparison to a method developed Each ITK transform now has either global support (affine trans- for ITK , via a standard brain labeling task. Finally, we discuss form) or local (or compact) support (a displacement field trans- future developments in the framework. form). If any map in a composite transform has global support then the composite transform has global support. Both “fixed” 2. MATERIALS AND METHODS and “moving” images may have initial transforms. This allows one 2.1. OVERVIEW OF THE UNIFIED FRAMEWORK to reduce “registration bias” that may be induced by asymmetric interpolation (Yushkevich et al., 2010). The overall purpose of the registration refactoring for ITK was to simplify the user experience and to accelerate and improve per- 2.1.4. Registration mappings should be applicable to a number of formance. Here, we summarize how ITK works toward these goals. popular data types, including DTI Our revisions to the ITK transform hierarchy validated and extended the ITK transforms for thread safety and applicability 2.1.1. Core Software Components to not only vectors but also tensors. Reorientation steps necessary Figure 1 sketches the ITK architecture at a high level. for diffusion tensor mappings are now included in ITK . Registration applications are known as “registration methods” as they were in ITK . The methods, with source contained in ITK ’s RegistrationMethodsv4 directory, hold together 2 4 All ITK metrics are set to be minimized. For instance, the the different subcomponents that make a working instantia- itkMattesMutualInformationImageToImageMetricv4 returns tion of a registration strategy. These subcomponents include the negative mutual information. Frontiers in Neuroinformatics www.frontiersin.org April 2014 | Volume 8 | Article 44 | 2 Avants et al. ITKv4 image registration ITKv4 Registration Framework FIGURE 1 | A schematic overview of the prototypical ITK registration metrics, greatly increasing efficiency; (4) automated (usually hidden) parameter method. This design is overall similar to that of ITK . A few key components estimators are available; (5) transforms may include high-dimensional differ: (1) optimizers require that transforms update themselves; (2) metrics and deformation fields. One additional difference (not shown) is that “fixed” images optimizers are multi-threaded; (3) memory is shared across both optimizers and may also have a transformation, although this is not modified by the optimizer. 2.1.5. Affine and deformable similarity metrics should look as (itkQuasiNewtonOptimizerv4), multiple objective opti- similar as possible mization (itkMultiGradientOptimizerv4), or global optimization (itkMultiStartOptimizerv4). The Metricsv4 frameworksupportsthisgoalinthatitisas trivial to implement a mutual information Demons algorithm as 2.1.8. GPU and multi-core acceleration will open up new it is to implement a sum of squared differences BSpline or affine applications for image registration registration algorithm. Thus, full plug-and-play support exists See GPUPDEDeformable for a GPU example. Furthermore, across transforms. the new metric framework N cores to accelerate metric, gradi- 2.1.6. Users should be able to combine multiple similarity metrics, ent and optimization steps. A recent real-world application of the some of which may operate on different data types new Insight ToolKit implementation of the symmetric normaliza- This is achievable with the existing itkMultiGradientOpt tion algorithm showed a speed-up of almost a factor of six when imizerv4 through the multivariate itkObjectToObjec comparing single core to eight core execution time. This speed-up tMultiMetricv4 or through the multi-channel traits (itkV is achieved by multi-threading the similarity metric, the gradi- ent descent update, the regularization terms and the composition ectorImageToImageMetricTraitsv4)that allow met- rics to deal with multi-channel pixels, all of which were con- components of the method. Thus, every essential step exploits intrinsic parallelism in the algorithm. Decreased execution time tributed for ITK .The itkObjectToObjectMultiMetri cv4 was used in our winning entry of the SATA 2013 “dog leg” means more rapid turnaround for users, faster turn-around in testing and higher throughput on large-scale computing tasks. challenge. 2.1.7. Optimizers and transformations should interact flexibly 2.1.9. Improve memory efficiency in optimization framework Optimizersv4 includes optimizers that are applicable to Memory optimizations are critical for efficient use of large local both linear and deformable transformations, which include transforms. In ITK , transform parameters are no longer copied convergence monitoring and enable 2nd order optimization within the optimizer, but rather left in-place in transform. Metric Frontiers in Neuroinformatics www.frontiersin.org April 2014 | Volume 8 | Article 44 | 3 Avants et al. ITKv4 image registration FIGURE 2 | Clockwise: Define x in  and z in  as the same material Second, we remove the shape (diffeomorphic) differences (D). Consequently, I J point but existing in different domains. The point y is in a domain that is we have a composite mapping, computed via the mutual information intermediate between  and  . The standard approach in the ITKv4 similarity metric, that identifies I(x) ≈ J(A(φ(x))) = J (y) = J(z). The I J mi Affine registration framework is to map image J (B) to image I (A) by first image J (y) represents J after application of the affine transformation A Affine identifying the linear transformation, →, between the images, shown in (C). i.e., J(A(x)). Code and data for this example are here. gradient memory is shared between optimizer and metric, and begin nomenclature definition modifications by the optimizer are done in place when possible. 4 A physical point: x ∈  where  is the domain, usually of Finally, we summarize ITK changes through quantitative an image. metrics: d n An image: I :  → R where n is the number of components per pixel and d is • Over 12 new multi-threaded image registration metrics are dimensionality. A second image is J. availableinv4. Domain map: φ :  →  where → may be replaced I J with any mapping symbol below. • Four application-level registration methods, with plug-and- Affine mapping: ↔ a low-dimensional invertible play architecture, are available for high-level inclusion in transformation: affine, rigid, translation, projects such as Slicer and SimpleITK. etc. • All contributions are unit-tested and have greater than 85% Affine mapping: → designates the direction an affine code coverage, in accordance with ITK standards. mapping is applied. • A complete refactoring of the ITK transform hierarchy that Deformation field:  deformation field mapping J to I.May makes transforms thread-safe, applicable to high-dimensional not be invertible. optimization and easily used in multi-core computing. Fourty- Spline-based mapping: e.g., B-Spline field mapping J to I. one classes, in total, were impacted by this refactoring. Diffeomorphism: Represented as, these are • We added transparent vector support to two key interpola- differentiable maps with differentiable tors that are used pervasively in ITK: the nearest neighbor and inverse. Ideally, the algorithm should output the inverse and forward mapping. linear interpolators. We added two new Gaussian interpolators. Composite mapping: φ = φ (φ (x)) is defined by→ where • An example of vector support for image metrics is in 1 2 φ is of type. itkMeanSquaresImageToImageMetricv4Vecto Not invertible:  indicates a mapping that is not rRegistrationTest.cxx. invertible. Perform image warping: As an example, → J represents the Below we will discuss: (0) an organizing nomenclature matched application of an affine transform → to to the ITK framework, (1) gradient-based optimization image J such that → J = J(A(x)). within the framework, (2) techniques to estimate optimiza- Similarity measure: or ≈ indicates the metric s that s s compares a pair of images. tion parameters for arbitrary metric and transformation end nomenclature definition combinations, (3) a ITK instance implementing general- ized diffeomorphic matching, (4) several applications of the representing algorithms and applications of registration. Ideally, updated framework within different neuroinformatics-relevant any standard algorithm can be written in the nomenclature below. domains. We would then write a standard Demons registration appli- cation that maps one image, J, into the space of I (presumably 2.2. NOMENCLATURE a template) as: The nomenclature below designates an image registration algorithm symbolically. This nomenclature is intended to be a descriptive and technically consistent system for visually I → J which symbolizes I ≈ J (A (φ (x))) , Frontiers in Neuroinformatics www.frontiersin.org April 2014 | Volume 8 | Article 44 | 4 Avants et al. ITKv4 image registration with A an affine mapping and φ a generic deformation. The nota- correlation and mutual information (Hermosillo et al., 2002). tion means that the algorithm first optimizes an affine mapping, Both are implemented in ITK for transformations with local ∂J(φ(x,p)) →,between J and I. This is followed by a deformation in the sec- and global support. The term is the gradient of J at ∂φ ond stage,,from → J to I. In terms of transformation compo- ∂φ φ(x)and is the Jacobian of the transformation taken with ∂p sition, we would write → J = J (x) = J(φ (φ (x))) w Affine Demons respect to its parameter. The transform φ(x, p) may be an affine where J is the result of warping J to I.The φ are the specific map i.e., φ(x, p) = Ax + t where A is a matrix and t a transla- functions corresponding to the schematic arrows. Note, also, that tion. Alternatively, it may be a displacement field where φ(x, p) = the tail of the arrow indicates the transform’s domain. The arrow- x + u(x)and u is a vector field. In ITK , both types of maps head indicates its range. Finally, we denote the similarity metric are interchangeable and may be used in a composite transform to as ≈ which indicates a sum of squared differences (the default compute registrations that map to a template via a schematic such similarity metric). ITK supports metrics such as mutual infor- ≈ ≈ as I ≈→ J, I → J, I  → J or, mixing similarity metrics, ≈ ≈ mi b cc mation, , or cross-correlation, . We will use this nomenclature mi cc I ≈ ≈ → J . cc mi i to write schematics for registration applications in the following The most commonly used optimization algorithm for image sections. registration is gradient descent, or some variant. In the above framework, the gradient descent takes on the form of 2.3. OPTIMIZATION FRAMEWORK The general ITK optimization criterion is summarized as: ∂M ∂M φ p , x = φ p + λ , ··· , , x , new old ∂p ∂p 1 n Find mapping φ(x, p) ∈ T such that M I, J,φ(x, p) is minimized. (1) where λ is the overall learning rate and the brackets hold the vector of parameter updates. While, for functional mappings, this formulation is not strictly In addition to basic gradient descent, we implement non- linear gradient descent optimization strategies which combine the correct, the practical implementation of even high-dimensional continuous transformations involves parameterization. The space conjugate gradient or gradient descent method with line search. In ITK , we implement the classic golden section approach to T restricts the possible transformations over which to optimize the mapping φ. The arguments to φ are its parameters, p,and identifying the optimal gradient step-size at each iteration. The generic conjugate gradient approach is performed via: the spatial position, x.Notethat, in ITK , the image I may also contain a mapping, although it is not directly optimized in most 2 2 cases. As will be seen later in the document, this mapping may ∇M −∇M t t − 1 γ = , (3) also be used within large deformation metrics. ∇M t − 1 The similarity metric, M, is perhaps the most critical com- CG =∇M + γ CG , t t t − 1 ponent in image registration. Denote a parameter set as p = (p , p ... p ). The metric (or comparison function between 1 2 n where CG is the conjugate gradient. The golden section line search images) is then defined by M(I, J,φ(x, p)). For instance, M = determines the specific weight,  , of the update to the current opt I(x) − J(φ(x, p)) i.e., the sum of squared differences (SSD) parameters such that metric. Its gradient with respect to parameter p is (using the chain rule), p = p + CG  . new old t opt ∂M ∂M ∂J φ(x, p) ∂φ Note that a naive application of gradient descent will not produce M = = | . (2) p x ∂p ∂J ∂φ ∂p i i a smooth change of parameters for transformations with mixed parameter types. For instance, a change, , to parameter p will This equation provides the metric gradient specified for sum of produce a different magnitude of impact on φ if p is a translation squared differences (at point x) but similar forms arise for the rather than a rotation. Thus, we develop an estimation framework that sets “parameter scales” (in ITK parlance) which, essentially, In the Demons example above, the reader may ask: why does the affine customize the learning rate for each parameter. The update to φ mapping, A, not appear “inside” the deformable mapping, φ? Indeed, this via its gradient may also include other steps (such as Gaussian ordering of transformations is feasible. However, this is not what we typically smoothing) that project the updated transform back to space T . use in our own practice of image registration and is not what we recommend. Multi-threading is achieved in the gradient computation, trans- The reason is that we usually seek a deformable mapping into a template formation update step and (if used) the regularization by dividing space that is common across many subjects (i.e., the “tail” is in the same the parameter set into computational units that correspond to domain across subjects). This enables methods such as Jacobian-based mor- phometry and other groupwise comparison conveniences. For example, see contiguous sub-regions of the image domain. dφ Figures 1 through 4 in (Kim et al., 2008) which explains the classical approach In termsofcode, theJacobian, | , is calculated at a dp of Jacobian-based morphometry as applied to traumatic brain injury. See physical point using the function ComputeJacobianWithRespect also (Gaser et al., 2001; Riddle et al., 2004; Dubb et al., 2005; Lepore et al., ToParameters(mappedFixedPoint, Jacobian).Notethatitiseval- 2006; Rohlfing et al., 2006; Studholme and Cardenas, 2007). Additional uses uated at point x notatpoint φ(x, p). We then use the func- of Jacobian-based morphometry are shown in our examples section in the “C” example and the “asymmetry” example. tion ComputeMovingImageGradientAtPoint(mappedMovingPoint, Frontiers in Neuroinformatics www.frontiersin.org April 2014 | Volume 8 | Article 44 | 5 Avants et al. ITKv4 image registration mappedMovingImageGradient) to compute the moving image which is regularized by the linear operator L. ITK uses Gaussian gradient when there is no pre-warping. ComputeMovingImageGra smoothing which is the Green’s kernel for generalized Tikhonov dientAtPoint uses central differences (or a gradient filter) in the regularization (Nielsen et al., 1997). This objective is readily opti- dJ(φ(x,p)) mized using an approach that is similar to that proposed by Beg. moving image space to compute the image gradient, . dφ If one is doing pre-warping, then we have an index access to Generalization of the LDDMM gradient for other metrics basi- cally follows (Hermosillo et al., 2002) with a few adjustments the warped moving image. We compute the warped image J as J (x) = J(φ(x, p)). Then, to accomodate diffeomorphic mapping. Figure 3 shows an ITK result on a standard example for large deformation registra- dJ dJ φ x, p d φ(x, p) tion. We will evaluate this diffeomorphic mapping, along with = (4) parameter estimation, in a later section. dx dφ dx −1 dJ φ(x, p) dJ d φ(x, p) 2.5. PARAMETER SCALE ESTIMATION = . We choose to estimate parameter scales by analyzing the result dφ dx dx of a small parameter update on the change in the magnitude of In code, we use ComputeMovingImageGradientAtIndex(index, physical space deformation induced by the transformation. The dJ mappedMovingImageGradient) to get and transform this impact from a unit change of parameter p may be defined in dx image gradient via the inverse Jacobian by calling mappedMovingI multiple ways, such as the maximum shift of voxels or the aver- mageGradient= TransformCovariantVector(mappedMovingImag age norm of transform Jacobians (Jenkinson and Smith, 2001). eGradient, mappedMovingPoint). Denote the unscaled gradient descent update to p as p.The goal is to rescale p to q = s · p,where s is a diagonal matrix 2.4. DIFFEOMORPHIC MAPPING WITH ARBITRARY METRICS diag(s , s ... s ), such that a unit change of q will have the same 1 2 n i The framework proposed above, in general form, encom- impact on deformation for each parameter i = 1 ... n. passes both classic affine mapping as well as more recent large As an example, we want φ(x, p ) − φ(x, p )= constant new old deformation strategies. Beg proposed the Large Deformation regardless of which of the i parameters is updated by the Diffeomorphic Metric Mapping (LDDMM) algorithm (Miller unit change. The unit is an epsilon value e.g., 1.e-3. Rewrite et al., 2005) which minimizes the sum of squared differences ∂J(φ(x,p)) ∂φ ∂φ ∂M ∂M ∂M , ··· , as , ··· , . To determine the ∂p ∂p ∂J ∂φ ∂p ∂p 1 n 1 n criterion between two images. LDDMM parameterizes a diffeo- relative scale effects of each parameter, p ,wecan factor out morphism through a time varying velocity field that is integrated 4 the constant terms on the outside of the bracket. Then the through an ODE. In ITK , we implement an alternative to ∂φ modified gradient descent step becomes diag(s) . We iden- ∂p LDDMM that also uses a time varying field and an ODE but tify the values of diag(s) by explicitly computing the values of minimizes a more general objective function: φ(x, p ) − φ(x, p ) with respect to an  change. A criti- new old cal variable, practically, is which x to choose for evaluation of φ(x, p ) − φ(x, p ). The corners of the image domain work new old E(v) = M I, J,φ + w Lv  dt. (5) 1,0 t well for affine transformations. In contrast, local regions of small radius (approximately 5) work well for transformations with local This is an instance of Equation (1) where w is a scalar weight and support. Additional work is needed to verify optimal parameters φ is a standard integration of the time-varying velocity field, v , for this new ITK feature. However, a preliminary evaluation is 1,0 t FIGURE 3 | An ITK diffeomorphic mapping of the type I  J. The “C” 15 5 integration steps were used in the Runge-Kutta ODE integration over and 1/2 “C” example illustrate the large deformations that may be achieved the velocity field. A two core MacBook Air computed this registration in 110 s. with time varying velocity fields. In this case, the moving (deforming) image The images each were of size 150 × 150. See C for a reproducible example of is the 1/2 “C.” The right panels illustrate the deformed grid for the this registration and the data. In addition, we provide an example of how the transformation of the “C” to 1/2 “C” (middle right) and its inverse mapping Jacobian determinant is computed from the deformation field resulting from (far right) which takes the 1/2 “C” to the reference space. The unit time this registration via an ANTs program CreateJacobianDeterminantIm interval is discretized into 15 segments in order to compute this mapping. age. Frontiers in Neuroinformatics www.frontiersin.org April 2014 | Volume 8 | Article 44 | 6 Avants et al. ITKv4 image registration performed in the results section. The new parameter scale esti- The ANTs contribution is valuable, in part, because of mation effectively reduces the number of parameters that the user the tremendous range of registration problems that exist in must tune from k + 1(λ plus the scales for each parameter type neuroinformatics and biomedical imaging in general. While it where there are k types) to only 1, the learning rate. is not possible to solve all registration problems with a gen- The learning rate, itself, may not be intuitive for a user to eral framework, one cannot afford to invent new solutions for set. The difficulty—across problem sets—is that a good learn- every instance one encounters. Our general optimization-driven strategies have proven to be invaluable to setting performance ing rate for one problem may result in a different amount of change per iteration in another problem. Furthermore, the dis- standards in a variety of application domains. In this section, crete image gradient may become invalid beyond one voxel. Thus, we highlight some of the lesser known capabilities of ANTs and it is good practice to limit a deformation step to one voxel spac- ITK with reproducible examples that include data and specific ing (Jenkinson and Smith, 2001). We therefore provide the users commands to ANTs and/or ITK. A list of these examples follows: the ability to specify the learning rate in terms of the maxi- mum physical space change per iteration. As with the parameter 1. The Basic Brain Mapping example shows how one can map scale estimation, the domain over which this maximum change two “whole head” T1-weighted MRIs where one is a tem- is estimated impacts the outcome and similar practices are rec- plate that contains a researcher’s prior knowledge defining ommended for both cases. This feature is especially useful for the “interesting” parts of the image. Within ITK ,this allowing one to tune gradient descent parameters without being domain knowledge is used to focus the Metricsv4 on concerned about which similarity metric is being used. That is, only those parts of the image while masking the remain- it effectively rescales the term λ∂M/∂p to have a consistent effect, der. Furthermore, a second part of this example shows how for a given λ,regardlessofthe metric choice.Incombination with the ITK composite transform may be used to initialize new our non-linear conjugate gradient approach (our current opti- registration solutions as well as how masking functionality mization of choice for linear registration), this strategy drastically may be employed to ignore information that is irrelevant or reduces the parameter setting burden for users. obstructive to the registration optimization. We have previ- ously employed these strategies in brain mapping with lesions 3. RESULTS (Avants et al., 2008; Kim et al., 2008, 2010, 2013; Tustison et 3.1. EXAMPLE APPLICATIONS OF THE ITK FRAMEWORK al., 2011). As part of our work in ITK refactoring, we built, in paral- 2. We use updated ITK methods in template construction with lel to library programming, an application interface that allows a reproducible example based on face and brain data: ANTs high-level access to the deep layers of ITK registration. These cur- template construction. This work has been employed in rently exist in the Advanced Normalization Tools (ANTs) software different species, age-ranges, and imaging modalities. The (link). While ANTs still serves as intermediate (vs. direct) access- resulting template is an image that captures the expected point to these tools, it provides a high-degree of customization shape and appearance as defined by the population sample, possibilities simply through a command line interface and script- transformation model and intensity comparison metric. ing. Therefore, a user is not required to write new low-level 3. A large deformation example implementing the classic “letter (interestingly, C++ is now considered “low-level”) software. C” example provided, originally, by Gary Christensen. While Despite its relative youth, the ANTs wrapping of ITK func- extremely flexible, these algorithms have not found a unique tionality has been employed with notable success in recent public, identity in terms of translational applications yet remain of unbiased, international evaluation studies. ANTs was instrumen- theoretical interest. This example shows a user how to define tal to a first-place finish in SATA 2013 in two of three categories the parameters of a registration based on optimizing a time (based on the median performance) where the ANTs approach varying velocity field. was considerably simpler than that employed by close finishers. 4. We present a separate example of how to compute landmark- While the evaluation of deep-gray matter registration showed based registration error. ITK uses LPS coordinates to rep- relatively subtle differences, the ANTs solution to the multivari- resent physical space. If you need to convert landmarks to ate canine leg data outclassed all other entrants. Notably, the physical space, see the discussion here: LPS physical space. ANTs solution used a multiple metric approach that simultane- We have an example illustrating how to change point coor- ously compared two modalities during registration as in Avants dinates and apply ITK transforms to landmarks here. This et al. (2008). In the cardiac data, the ANTs solution was the exercise can be useful for landmark-based registration or in only one that was fully automated resulting in a ≈15% per- evaluating registration accuracy. formance loss which can easily be overcome by a modicum of 5. We show how to perform motion correction to time series user intervention. Furthermore, ANTs/ITK-based methods fin- data here although we do not claim this approach is opti- ished a clear first-place in the BRATS 2013 challenge. Our entry mal. The method registers each frame from a 4D time series used intensity asymmetry as a key feature to segment brain to a fixed reference image and stores the resampled set in a tumors based on multiple modality MRI. Thus, these methods new 4D image. All transformation parameters are stored in a are within the leading ranks of image registration methodolo- corresponding csv file. gies as evaluated in recent work as well as in the more traditional 6. An advanced example with heavy use of statistics via R and brain (Klein et al., 2010) and lung CT (Murphy et al., 2011) ANTsR is in a study of public test-retest fmri data. This study studies. is not published and may be subject to change. Frontiers in Neuroinformatics www.frontiersin.org April 2014 | Volume 8 | Article 44 | 7 Avants et al. ITKv4 image registration 7. The classic car example shown in ANTs talks is here. This 3.2.1. Parameter estimation across metrics illustrates the benefit of mutual information in deformably ITK provides similarity metrics that may be applied for both mapping wild-type images and highlights the fact that deformable and affine registration. In a previous section, we pro- ITK applications exist outside of medical imaging. vided a parameter estimation strategy that is applicable to both 8. A basic multistart example. A more advanced example deformable and affine transformations with arbitrary metrics. for brain mapping with a template mask is available in Denote images I, J, K, where the latter two are “moving” images, antsCorticalThickness.sh. These optimization methods over- and K is an intensity-inverted version of J.Wethenevaluatethe come local minima by running registration searches from following schema, a variety of starting points and greedily storing the best solution. I ≈→ J, I ≈ → K, I ≈ → K cc mi 9. This asymmetry analysis example uses the “virtual domain” where, for each schematic, we use the corresponding metric for feature to reduce bias caused by mapping an image asymmet- both affine and diffeomorphic mapping. Furthermore, we keep rically to a reference. Note that we measure the point-wise the same parameters for each registration by exploiting parameter asymmetry, in this example, via the Jacobian determinant scale estimators. Figure 5 shows the candidate images for this test. image as in Kim et al. (2008). If one repeats this analy- As shown in Figure 5, very similar results are achieved for each sis across a population—and maps the Jacobian measure- schematic without additional parameter tuning. To determine ments of asymmetry to a common space—then one may this quantitatively, we perform registration for each schematic perform a statistical analysis of population-level asymme- and then compare the Dice overlap of a ground-truth three-tissue try. Longitudinal analysis and asymmetry analysis potentially segmentation. For each result, we have the Dice overlap of dark suffer from the same confound Yushkevich et al. (2010). A tissue (cerebrospinal fluid, CSF), medium intensity tissue (gray related longitudinal mapping script is here. matter) and bright tissue (white matter). For the mean squares 10. We also show how manual labelings can be used to restrict metric, we have: 0.588, 0.816, and 0.90; for CC, we have: 0.624, registration in a challenging registration scenario (this exam- 0.786, 0.882; for MI, we have: 0.645, 0.779, 0.858. Mutual infor- ple will be improved in the future): registration guided by mation does best for the CSF while mean squares does best for (crude) labels. other tissues. CC performs in the mid-range for all classes of 11. A simple orange to apple RGB image registration example for tissue. Thus, a single set of tuned parameters provides a reason- color images is listed at: itkMeanSquaresImageToImageMet- able result for an affine plus diffeomorphic mapping across three ricv4VectorRegistrationTest. If one compiles the ITK tests, different metrics. While improvement might be gained by fur- then this example can be run to produce Figure 4. ther tuning for each metric, this result shows that our parameter estimation method achieves the goal of reducing user burden. These examples cover many applications for which no “ground 3.2.2. Automated Brain Labeling Task truth” evaluation data exists. The next section seeks to add some quantitative reference to these examples. First, we show flexibility All R and bash analysis scripts for this section are here: https://github.com/stnava/ITKv4Documentation/tree/frontiers/sc and consistency of our framework in a simple example comparing registration with a variety of metrics and a consistent parameter ripts. The ITK core functionality formed the heart of the ref- erence results provided for the SATA2013 challenge at MICCAI set. Second, we quantify the benefit of ITK registration in com- parison to a method implemented based upon ITK registration 2013. In this sense, these methods have been heavily evaluated on both basic brain mapping challenges (SATA2013’s dien- technology. cephalon challenge in which ITK -based methods finished first), multivariate registration challenges (the canine MRI / dog 3.2. EVALUATION leg challenge of SATA2013 in which ITK -based methods were We first investigate the ability of our automated parameter esti- overwhelmingly the top finisher) and in the cardiac challenge mation to facilitate parameter tuning across metrics. Second, we 4 3 (in which ITK -based methods were the only fully automated compare ITK and ITK registration implementations with approach). However, for completeness, we provide an additional respect to a standard automated brain labeling task. evaluation here which focuses on comparison to a ITK method, BRAINSFit, in a different dataset than previously used to evaluate ANTs or ITK . As ground truth, we use T1 MRI data from 33 2-year old subjects as described in Gousias et al. (2008) and available at http://www.brain-development.org. Each subject’s brain is man- ually parcellated into 83 distinct regions that include ventricles, cortical areas, white matter and deep gray matter regions such as the amygdala, hippocampus and thalamus. One benefit of this data is that some of these anatomical regions are rela- FIGURE 4 | This RGB image registration example employs ITK code tively easy to align (the caudate) whereas others are relatively that repurposes a scalar metric (itkMeanSquaresImageToImageMetr difficult to align due to their small size (amygdala) or incon- icv4) for multichannel registration. sistent shape across subjects (the inferior frontal gyrus). Thus, Frontiers in Neuroinformatics www.frontiersin.org April 2014 | Volume 8 | Article 44 | 8 Avants et al. ITKv4 image registration FIGURE 5 | Three reference images, I (left), J (middle top), and K (right affine and deformable case. Thus, two parameters had to be optimized. We top), are used to illustrate the robustness of our parameter scale then applied these same parameters to register I and K via both correlation estimation for setting consistent parameters across both metrics and and mutual information. The resulting registrations (bottom row) were all of transform types. K is the negation of J andisusedtotestthe correlation similar quality. Further, the same metric is used for both affine and and mutual information registrations. We optimized, by hand, the step-length diffeomorphic mapping by exploiting the general optimization process given parameters for one metric (the sum of squared differences) for both the in Equation (1). FIGURE 6 | We compare a ITK composite schema as I ≈ ≈ → J example ANTs-based large-deformation result from the dataset is shown for cc mi i for mapping a set of {J } images to a template I to a ITK schema: illustration where we render the extracted brains as well as show select axial I ≈  ≈ → J . We use this schematic in a registration-based slices. All registrations were run on the original MRI data with no mi b mi i segmentation of multiple brain structures in a pediatric population as a preprocessing except what is done by ANTs or BRAINSFit internally. Overlap benchmark for algorithm performance, similar to Klein et al. (2010). An improvement from v3 to v4, quantified via paired t-test, is highly significant. Frontiers in Neuroinformatics www.frontiersin.org April 2014 | Volume 8 | Article 44 | 9 Avants et al. ITKv4 image registration we anticipate that performance gains due to new technology in contained within the ANTs repository. The git hashtag for the ITK will be most prominent in the more variable and chal- BRAINSFit version https://github.com/BRAINSia/BRAINSTools lenging regions. Figure 6 summarizes the study nomenclature used in this evaluation is ad7e114ab1c92bd800819b80e054825 and shows a single image pair selected from this data along 9398931c8. Both programs were run on a MacBook Pro running with the registration result given by ITK . Figure 7 summa- OS X 10.9 (13A3028) with a 2.6 GHz Intel Core i7 and 16 GB of rizes these evaluation results. The scripts for running this study RAM. are available at https://github.com/stnava/ITKv4Documentation/ To help isolate which subject pairs to deformably register, tree/frontiers/scripts. The git hashtag for the ANTs version used we first clustered the initial dataset based on an all pairs affine in this evaluation https://github.com/stnava/ANTs is ce8b5a741 registration which revealed five representative subject clusters. 4ae9e389071d756c5f36ee6cecbcfd8. The associated ITK tag is The subject pairings used during evaluation were chosen such FIGURE 7 | Above, a barplot shows the mean Dice score for each region ITK SyN algorithm, with its classic neighborhood correlation metric, and each algorithm, sorted by ANTs performance. Below, we use star outperforms BRAINSFit in several regions and more strongly in some subject plots of per-brain-region Dice overlap to compare, for each subject, the pairs than others. The legend for the plots is at lower right and shows the 4 3 ITK implementation of SyN with the ITK -based BRAINSFit algorithm. The maximum possible value for each region. Frontiers in Neuroinformatics www.frontiersin.org April 2014 | Volume 8 | Article 44 | 10 Avants et al. ITKv4 image registration Table 1 | Dice overlap for ANTs and BRAINSFit where only regions that each subject pair contained the most representative subject with q-value < 0.01 are shown. for one cluster paired with the most representative subject from another cluster; thus, the study design allows us to focus, in a meanants meanbfit FDR-adjusted principled manner, on a set of representative shape comparisons p-value where the comparisons are made across different image types. The new methods in ITK show enhanced performance within R Hippocampus 0.835 0.770 0.0083 all registration pairs. The mean overall performance gain was L Hippocampus 0.834 0.782 0.0063 approximately 6.3% with a standard deviation of 5% and T- R Amygdala 0.843 0.752 0.0027 statistic/p-value, over all structures, of 12.6 / p < 1.e − 16. We L Amygdala 0.851 0.773 0.0037 also identified which regions were most improved in ITK vs. L AnteriorTemporalLobeMe 0.849 0.752 0.0005 ITK . These regions include the left and right insula, the brain- dialPart stem, the superior temporal gyrus, parahippocampal gyrus, puta- R AnteriorTemporalLobeLate 0.754 0.633 0.0002 men, and the substantia nigra. Table 1 lists all structures and ralPart the mean Dice score for each algorithm, along with the p-value. L AnteriorTemporalLobeLat 0.786 0.668 0.0083 Figure 7 summarizes all of these findings by using star plots to eralPart R Gyri parahippocampalis et 0.813 0.756 0.0037 visualize the Dice results for every region in every subject. ambiens We also recorded the amount of processing time spent on each L Gyri parahippocampalis et 0.823 0.749 0.0007 subject, for each algorithm. Noting that the ITK algorithm also ambiens provides a dense and high-resolution forward and inverse trans- R Superior temporal gyrus 0.892 0.827 0.0001 form and does explicit transformation regularization to guarantee posterior a diffeomorphism, the algorithm takes, on average, 5 times as long 3 R Medial and inferior 0.891 0.823 0.0007 as the ITK BRAINSFit algorithm (≈10 min), assuming default temporal gyri settings. Much of this time, in ANTs, is taken up by full resolution R Lateral occipitotemporal 0.801 0.690 0.0023 image registration. If this fine level is avoided, then the disparity gyrus (gyrus fusiformis) in timing reduces to less than a factor of two, without much loss R Cerebellum 0.953 0.937 0.0066 in accuracy. Note also that ANTs and BRAINSFit each use a differ- Brainstem 0.944 0.923 0.0001 ent multithreading strategies, similarity metric implementations, R Insula 0.910 0.864 0.0001 rigid/affine registration mechanisms and optimizers making this L Insula 0.915 0.863 0.0001 overall comparison less than ideal. L Occipital lobe 0.895 0.819 0.0049 L Posterior temporal lobe 0.912 0.862 0.0001 4. DISCUSSION L Parietal lobe 0.885 0.815 0.0052 ITK is a community built and maintained toolkit and is a public R Putamen 0.909 0.869 0.0002 resource for reproducible methods. The updated ITK registra- L Putamen 0.911 0.862 0.0003 tion framework provides a novel set of user-friendly parameter R Pallidum 0.838 0.721 0.0037 setting tools and benchmark implementations of both standard L Pallidum 0.846 0.766 0.0049 and advanced algorithms. Robustness with respect to parameter L Lingual gyrus 0.876 0.795 0.0037 settings has long been a goal of image registration and ITK takes L Cuneus 0.825 0.684 0.0083 valuable steps toward the direction of automated parameter L Lateral orbital gyrus 0.728 0.598 0.0037 selection. The primary decision left up to the user, currently, L Substantial nigra 0.831 0.741 0.0023 is the feature scale at which registration should be performed. L Superior temporal gyrus 0.818 0.752 0.0043 E.g., whether the registration should focus on coarse features, anterior part fine features, etc and the different resolutions at which this should be done. While we have provided a reproducible reference comparison of registration-based brain labeling in this paper, we implementations for speed and memory. In time, it may be intend to have a more extensive series of benchmark performance possible to extend the design philosophy used here to GPU studies completed on datasets beyond the brain. However, the implementations. However, our ability to interface low and number of possible applications exceeds what can possibly be high-dimensional transformations depends heavily on generic evaluated by the core ITK developer community. Community programming. This style is less well-developed (and less well involvement is needed in order to increase the number of possible understood) in GPU applications which depend, to some registration applications and metric/transform/optimizer/data extent, on specialization. The current framework is amenable combinations that have been evaluated. At the same time, to groupwise registration strategies when used in combina- documentation, usability and examples must be pro- tion with a computing cluster. However, single core group- vided by the development team in order to improve user wise strategies are not currently implemented although one involvement. may consider basing an implementation on exisiting multi- 4.1. FUTURE WORK metric/multivariate registration tools within the current code Future work will enhance the depth and breadth of docu- base. While ITK does contain a statistics infrastructure, we mentation as well as seek to further optimize the current currently prefer using R and ANTsR for analyzing our data. Frontiers in Neuroinformatics www.frontiersin.org April 2014 | Volume 8 | Article 44 | 11 Avants et al. ITKv4 image registration However, the lack of visualization methods in ITK means volumetry of brain lateral ventricles in mri. Neuroimage 13(6 pt 1), 1140–1145. doi: 10.1006/nimg.2001.0771 that one must still move to another package to look at Gousias, I. S., Rueckert, D., Heckemann, R. A., Dyet, L. E., Boardman, one’s results. Therefore, direct interfaces to R remain useful. J. P., Edwards, A. D., et al. (2008). Automatic segmentation of brain mris SimpleITK also has a promising R interface that is similar to of 2-year-olds into 83 regions of interest. Neuroimage 40, 672–684. doi: ANTsR. 10.1016/j.neuroimage.2007.11.034 A primary challenge to the future of ITK includes, beyond Hermosillo, G., Chefd’Hotel, C., and Faugeras, O. (2002). A variational approach to multi-modal image matching. Intl.J.Comp. Vis. 50, 329–343. doi: documentation, reduced C ++ fluency. As ITK leverages 10.1023/A:1020830525823 several advanced features of C ++, even experienced devel- Jack, C. R. Jr., Knopman, D. S., Jagust, W. J., Shaw, L. M., Aisen, P. S., opers may find it difficult to contribute meaningfully to the Weiner, M. W., et al. (2010). Hypothetical model of dynamic biomark- ITK software base. Therefore, the ITK community must also ers of the Alzheimer’s pathological cascade. Lancet Neurol. 9, 119–128. doi: 10.1016/S1474-4422(09)70299-6 seek to educate potential future contributors not only on ITK Jenkinson, M., and Smith, S. (2001). A global optimisation method for robust affine but also, at times, on the fundamentals or advanced exten- registration of brain images. Med. Image Anal. 5, 143–156. doi: 10.1016/S1361- sions of C ++. A second major hurdle is that ITK includes 8415(01)00036-6 a host of generic registration ingredients. However, many of Kikinis, R., and Pieper, S. (2011). 3d slicer as a tool for interactive brain tumor the most compelling new application domains require special- segmentation. Conf. Proc. IEEE Eng. Med. Biol. Soc. 2011, 6982–6984. doi: 10.1109/IEMBS.2011.6091765 ization. Specialization may be needed for a specific imaging Kim, J., Avants, B., Patel, S., Whyte, J., Coslett, B. H., Pluta, J., et al. (2008). modality, via hardware interface or in the use of domain- Structural consequences of diffuse traumatic brain injury: a large defor- specific prior knowledge. Therefore, we envision the next phase mation tensor-based morphometry study. Neuroimage 39, 1014-1026. doi: of ITK development may focus on using the toolkit to sup- 10.1016/j.neuroimage.2007.10.005 port its specialization in solving high-impact and translational Kim, J., Avants, B., Whyte, J., and Gee, J. C. (2013). Methodological considerations in longitudinal morphometry of traumatic brain injury. Front. Hum. Neurosci. applications. Hopefully, this transition will occur in the near 7:52. doi: 10.3389/fnhum.2013.00052 future. Kim, J., Whyte, J., Patel, S., Avants, B., Europa, E., Wang, J., et al. (2010). Resting cerebral blood flow alterations in chronic traumatic brain injury: an arte- rial spin labeling perfusion fmri study. J. Neurotrauma 27, 1399-1411. doi: REFERENCES 10.1089/neu.2009.1215 Ackerman, M. J., and Yoo, T. S. (2003). The visible human data sets (VHD) and insight toolkit (ITk): experiments in open source software. AMIA Annu. Symp. Klein, S., Staring, M., Murphy, K., Viergever, M. A., and Pluim, J. P. W. (2010). elastix: a toolbox for intensity-based medical image registra- Proc. 2003:773. Avants, B. B., Tustison, N. J., Song, G., Cook, P. A., Klein, A., and Gee, J. C. (2011). tion. IEEE Trans. Med. Imaging 29, 196–205. doi: 10.1109/TMI.2009. A reproducible evaluation of ANTs similarity metric performance in brain image registration. Neuroimage 54, 2033–2044. doi: 10.1016/j.neuroimage. Lepore, N., Brun, C. A., Chiang, M.-C., Chou, Y.-Y., Dutton, R. A., Hayashi, K. M., et al. (2006). Multivariate statistics of the jacobian matrices in tensor based mor- 2010.09.025 phometry and their application to hiv/aids. Med. Image Comput. Comput. Assist. Avants, B., Duda, J. T., Kim, J., Zhang, H., Pluta, J., Gee, J. C., et al. (2008). Multivariate analysis of structural and diffusion imaging in traumatic brain Interv. 9(pt 1), 191–198. Metz, C. T., Klein, S., Schaap, M., van Walsum, T., and Niessen, W. J. (2011). injury. Acad. Radiol. 15, 1360-1375. doi: 10.1016/j.acra.2008.07.007 Baloch, S., and Davatzikos, C. (2009). Morphological appearance mani- Nonrigid registration of dynamic medical imaging data using nd + t b-splines and a groupwise optimization approach. Med. Image Anal. 15, 238–249. doi: folds in computational anatomy: groupwise registration and morphologi- cal analysis. Neuroimage 45(1 Suppl), S73–S85. doi: 10.1016/j.neuroimage. 10.1016/j.media.2010.10.003 Miller, M. I., Beg, M. F., Ceritoglu, C., and Stark, C. (2005). Increasing the power 2008.10.048 Bearden, C. E., van Erp, T. G. M., Dutton, R. A., Tran, H., Zimmermann, L., Sun, of functional maps of the medial temporal lobe by using large deformation dif- feomorphic metric mapping. Proc. Natl. Acad. Sci. U.S.A. 102, 9685–9690. doi: D., et al. (2007). Mapping cortical thickness in children with 22q11.2 deletions. Cereb. Cortex 17, 1889–1898. doi: 10.1093/cercor/bhl097 10.1073/pnas.0503892102 Murphy, K., van Ginneken, B., Reinhardt, J. M., Kabus, S., Ding, K., Deng, X., Buckheit, J. B., and Donoho, D. L. (1995). “WaveLab and reproducible research,” in Wavelets and statistics,eds A. Antoniadis andG.Oppenheim (New York,NY: et al. (2011). Evaluation of registration methods on thoracic CT: the EMPIRE10 challenge. IEEE Trans. Med. Imaging 30, 1901–1920. doi: 10.1109/TMI.2011. Springer), 55–81. Chen, M., Lu, W., Chen, Q., Ruchala, K. J., and Olivera, G. H. (2008). A simple 2158349 fixed-point approach to invert a deformation field. Med. Phys. 35, 81–88. doi: Nielsen, M., Florack, L., and Deriche, R. (1997). Regularization, scale- space, and edge detection filters. J. Math. Imaging Vis. 7, 291–307. doi: 10.1118/1.2816107 Cheung, M. R., and Krishnan, K. (2009). Interactive deformation registration of 10.1023/A:1008282127190 Peyrat, J.-M., Delingette, H., Sermesant, M., Xu, C., and Ayache, N. (2010). endorectal prostate mri using itk thin plate splines. Acad. Radiol. 16, 351–357. doi: 10.1016/j.acra.2008.09.011 Registration of 4d cardiac ct sequences under trajectory constraints with multi- channel diffeomorphic demons. IEEE Trans. Med. Imaging 29, 1351–1368. doi: Christensen, G. E., Rabbitt, R. D., and Miller, M. I. (1996). Deformable templates using large deformation kinematics. IEEE Trans. Image Process. 5, 1435–1447. 10.1109/TMI.2009.2038908 Riddle, W. R., Li, R., Fitzpatrick, J. M., DonLevy, S. C., Dawant, B. M., doi: 10.1109/83.536892 Dubb, A., Xie, Z., Gur, R., Gur, R., and Gee, J. (2005). Characterization of brain and Price, R. R. (2004). Characterizing changes in mr images with color- coded jacobians. Magn. Reson. Imaging 22, 769–777. doi: 10.1016/j.mri. plasticity in schizophrenia using template deformation. Acad. Radiol. 12, 3–9. doi: 10.1016/j.acra.2004.06.009 2004.01.078 Rohlfing, T., Sullivan, E. V., and Pfefferbaum, A. (2006). Deformation-based Fedorov, A., Li, X., Pohl, K. M., Bouix, S., Styner, M., Addicott, M., et al. (2011). Atlas-guided segmentation of vervet monkey brain MRI. Open Neuroimag. J. 5, brain morphometry to track the course of alcoholism: differences between intra-subject and inter-subject analysis. Psychiatry Res. 146, 157–170. doi: 186–197. doi: 10.2174/1874440001105010186 10.1016/j.pscychresns.2005.12.002 Floca, R., and Dickhaus, H. (2007). A flexible registration and evalua- tion engine (f.r.e.e.). Comput. Methods Programs Biomed. 87, 81–92. doi: Rueckert, D., Sonoda, L. I., Hayes, C., Hill, D. L., Leach, M. O., and Hawkes, D. J. (1999). Nonrigid registration using free-form deformations: applica- 10.1016/j.cmpb.2007.04.009 Gaser, C., Nenadic, I., Buchsbaum, B. R., Hazlett, E. A., and Buchsbaum, M. S. tion to breast mr images. IEEE Trans. Med. Imaging 18, 712–721. doi: 10.1109/42.796284 (2001). Deformation-based morphometry and its relation to conventional Frontiers in Neuroinformatics www.frontiersin.org April 2014 | Volume 8 | Article 44 | 12 Avants et al. ITKv4 image registration Shelton, D., Stetten, G., Aylward, S., Ibáñez, L., Cois, A., and Stewart, Yushkevich, P. A., Avants, B. B., Das, S. R., Pluta, J., Altinay, M., Craige, C. (2005). Teaching medical image analysis with the insight C., et al. (2010). Bias in estimation of hippocampal atrophy using toolkit. Med. Image Anal. 9, 605–611. doi: 10.1016/j.media.2005. deformation-based morphometry arises from asymmetric global normal- 04.011 ization: an illustration in adni 3 t mri data. Neuroimage 50, 434-445. doi: Studholme, C., and Cardenas, V. (2007). Population based analysis of directional 10.1016/j.neuroimage.2009.12.007 information in serial deformation tensor morphometry. Med. Image Comput. Comput. Assist. Interv. 10(pt 2), 311–318. Conflict of Interest Statement: The authors declare that the research was con- Tustison, N., Avants, B., Cook, P., Kim, J., Whyte, J., Gee, J., et al. (2011). ducted in the absence of any commercial or financial relationships that could be “Multivariate analysis of diffusion tensor imaging and cortical thickness maps construed as a potential conflict of interest. in a traumatic brain injury (tbi) cohort using advanced normalization tools (ants) Journal of Neurotrauma 28 A111,” in 29th Annual National Neurotrauma Received: 01 November 2013; accepted: 30 March 2014; published online: 28 April Symposium (Hollywood Beach, FL). 2014. van Dalen, J. A., Vogel, W., Huisman, H. J., Oyen, W. J. G., Jager, G. J., Citation: Avants BB, Tustison NJ, Stauffer M, Song G, Wu B and Gee JC (2014) The and Karssemeijer, N. (2004). Accuracy of rigid CT-FDG-PET image registra- Insight ToolKit image registration framework. Front. Neuroinform. 8:44. doi: 10.3389/ tion of the liver. Phys. Med. Biol. 49, 5393–5405. doi: 10.1088/0031-9155/49/ fninf.2014.00044 23/014 This article was submitted to the journal Frontiers in Neuroinformatics. Vercauteren, T., Pennec, X., Perchant, A., and Ayache, N. (2009). Diffeomorphic Copyright © 2014 Avants, Tustison, Stauffer, Song, Wu and Gee. This is an open- demons: efficient non-parametric image registration. Neuroimage 45(1 Suppl), access article distributed under the terms of the Creative Commons Attribution S61–S72. doi: 10.1016/j.neuroimage.2008.10.040 License (CC BY). The use, distribution or reproduction in other forums is permit- Yoo, T. S., Ackerman, M. J., Lorensen, W. E., Schroeder, W., Chalana, V., Aylward, ted, provided the original author(s) or licensor are credited and that the original S., et al. (2002). Engineering and algorithm design for an image processing api: publication in this journal is cited, in accordance with accepted academic prac- a technical report on itk–the insight toolkit. Stud. Health Technol. Inform. 85, tice. No use, distribution or reproduction is permitted which does not comply with 586–592. these terms. Frontiers in Neuroinformatics www.frontiersin.org April 2014 | Volume 8 | Article 44 | 13

Journal

Frontiers in NeuroinformaticsPubmed Central

Published: Apr 28, 2014

There are no references for this article.