Neuroimage 50(4): 1485-1496

PMID: 20097291

Attribute Vector Guided Groupwise Registration

1. Introduction

Medical image registration has been a hot area of research for decades, due to its valuable applications to clinical studies for comparison of longitudinal and cross-sectional data, computer-assisted diagnosis and treatment guidance, and quantitative monitoring of the disease progression (Hill et al., 2001; Maintz and Viergever, 1998; Zitová and Flusser, 2003). For example, to analyze and monitor tumor growth based on the observation of patient images at different time points, meaningful diagnosis can only be conducted after aligning those images to a common space and preserving morphological variability related to the disease (Hajnal et al., 1995; Holden et al., 2000; Woods et al., 1992).

Pairwise image registration seeks to register the moving image to the space of the fixed image by minimizing the distance between two images. For example, given a set of corresponding landmarks on both the moving and the fixed images, the minimal landmark distance can be measured and optimized using the Iterative Closest Point (ICP) algorithm (Maintz and Viergever, 1998). However, no matter done manually or automatically, the determination of correspondences between these landmarks is challenging and usually ends in unreliable results. To this end, voxel-based registration methods are proposed in the literature, where energy terms capturing the distance between two intensity images are minimized in registration. There are a variety of image distance/similarity metrics available nowadays, ranging from the L^{norm of the intensity differences between two images (}Christensen et al., 1996), to information theoretic measurements such as mutual information (MI) (Maes et al., 1997; Wells et al., 1996) and normalized mutual information (NMI) (Studholme et al., 1999) which are both capable of measuring multi-modality similarity. On the other hand, transformation models also play critical roles in registration. Non-linear transformations, which have much higher degrees of freedom (DOF) than linear counterparts such as rigid or affine transformations, can generally better describe the subtle changes of anatomical structures in medical images. Such transformations can usually be represented by a linear combination of polynomials (Woods et al., 1998a; Woods et al., 1998b), basis functions (Friston et al., 1995), or B-Splines (Rueckert et al., 1999), etc.

To better understand and analyze the group similarity and variation within a population, it is now more important that a group of images, instead of a single pair, is processed and evaluated together in clinical and research applications. However, if pairwise registration methods are directly applied to more than two images, one image needs to be selected as the template and would inevitably bias the subsequent image analysis. To this end, several groupwise registration algorithms have been proposed recently to seek for spatial correspondences among a group of subjects and simultaneously warp them to a common space (Crum et al., 2004; Toga and Thompson, 2001). For example, Seghers et al. (Seghers et al., 2004) proposed a groupwise registration algorithm which performs exhaustive pairwise registrations between all possible pairs of images in the group. The final atlas is constructed by averaging all the images after they have been warped to their mean morphological templates, each of which is determined by averaging the transformations between the given image and all other images. In the method proposed by Park et al. (Park et al., 2005), the image closest to the population mean geometry, which is estimated through multi-dimensional scaling (MDS) of all images in the beginning, is defined as the tentative template. All images are registered to the template, as the template is refined by updating itself with the mean of those aligned images. Recently, minimum spanning tree (MST) was introduced to describe the linkages of different subjects within a group (Munsell et al., 2009), and the tree structure was used to reduce the groupwise registration into a series of pairwise registration problems. Though the method focuses on 3D shapes, it could easily be extended to medical images. While enjoying effectiveness of groupwise registration, the methods described above suffer from high computational complexity, especially when the size of the image population increases.

More efficient groupwise registration methods, which completely avoid the traditional pairwise registration paradigm, have consequently been investigated. In Joshi et al. (Joshi et al., 2004), a widely accepted large deformation diffeomorphic registration framework (Beg et al., 2005) is adapted for groupwise registration. In Zöllei et al. (Zöllei et al., 2005), a groupwise registration scheme based on the congealing framework (Learned-Miller, 2006; Miller et al., 2000) is proposed, where intensity based entropy drives a gradient-based stochastic optimizer and pushes all images to the population center simultaneously. This method, which works under the constraint of affine transformations only in Zöllei et al., is further extended in Balci et al. (Balci et al., 2007) by using B-Splines to model nonlinear deformations. Stack entropy, the cost function for the formulation in Balci et al. and Zöllei et al., measures the compactness of the cross-subject intensity distribution for a certain voxel location. If all images were perfectly registered, the stack entropy would ideally approach to zero. However, the shortcomings of these methods are that they rely on voxelwise intensity information, leaving out regional contextual information and anatomical shape characteristics, and that each voxel contributes equally to the registration.

To improve the robustness of registration, we propose a novel groupwise registration algorithm in this paper, which utilizes multi-scale attribute vectors, rather than only image intensities, as feature descriptors. To better capture the variability of each attribute in the attribute vector, a local neighborhood, instead of a single voxel, is taken into consideration. The neighborhood works as a local pattern, and is described by the distribution of attributes within the neighborhood. The inhomogeneity of corresponding local patterns across all input images is measured by the Jensen-Shannon (JS) divergence (Lin, 1991), whose overall sum for all attributes and for all voxels is used as the cost function for groupwise registration. We provide a gradient-based solution to optimize this cost function. In the end, all subjects in the group will be gradually registered to a common space without introducing any bias.

Our groupwise registration method was evaluated on both real data and simulated data. For the real data, we performed our method on the NIREP NA0 dataset of 16 brain images and compared the overlap ratios on 32 manually delineated ROIs. In another experiment, we simulated atrophies at pre-central gyrus (PCG) and superior temporal gyrus (STG). After normalizing both the atrophic subjects and the normal controls, the Jacobian determinant maps of the estimated deformation fields were used for t-test in SPM to measure the atrophy discrepancy capability. All experimental results were compared with the intensity guided groupwise registration method which was implemented in Balci et al. (Balci et al., 2007). And the comparison demonstrated that our method could achieve more accurate and consistent groupwise registration.

2. Methods

One novelty of our groupwise registration algorithm is that we utilize multi-scale attribute vectors as voxel anatomical signatures to guide groupwise registration. For each voxel location in the image domain, we consider the statistical information of attributes in its local neighborhood. Specifically, a set of probability density functions (PDFs) is estimated to describe the corresponding local patterns from individual images. We then use the JS divergence to measure inhomogeneity of the cross-subject PDFs. By minimizing the overall JS divergence integrated for all attributes and all voxels, a collection of transformations are optimized to pull each subject from its own space, which is initially far away from each other, to the final common space. We will describe our method in more details in the following sections.

2.1 Attribute and Attribute Vector

In medical image analysis, it is important to build anatomical models that take the underlying anatomy into account, rather than only the similarity on image intensities. This motivates the utilization of descriptive features which are able to more effectively represent different shapes of complex anatomical structures in medical images. This is in the spirit of HAMMER (Shen and Davatzikos, 2002), where geometric invariant moments are introduced as elements to form the attribute vector, and elastic registration is guided by those hierarchically selected driving voxels whose correspondences can be estimated more reliably based on their attribute vectors. These driving voxels are usually located at transition areas between two types of tissue labels (e.g., the interface of white matter and gray matter), and highlight the importance of contextual information to registration.

The composition of the attribute vector is quite flexible and can be easily fine-tuned to cater for different applications. Besides complex feature descriptors that are popular in the computer vision society, relatively simple attributes can also contribute to achieve high registration accuracy. Under the assumption that intensity transitions should occur at the same location after two images are well aligned, the metric of the registration between two images is then related to the directions of the intensity gradients in Haber and Modersitzki (Haber and Modersitzki, 2006). Recently, a combination of features, including the gradients and the Hessians in the intensity space, is applied to the registration of cervical MR images in Staring et al. (Staring et al., 2009). Although the studies mentioned above belong to the category of pairwise registration, they have indeed confirmed that the gradient information is quite helpful to guide registration. Therefore, better registration results can be expected with the help of more descriptive attributes.

For N images under groupwise registration, we denote the attribute vector at voxel location x of the i-th image as a D-tuple ${\overset{⇀}{A}}_{i, x} ≜ (a_{i, x}^{1}, \dots, a_{i, x}^{j}, \dots, a_{i, x}^{D})$ , where D denotes the number of attributes. Each attribute $a_{i, x}^{j}$ is a scalar value, and describes a type of features which can constitute not only the voxel intensity but also more regional information (as ① in Fig. 1). Therefore, the groupwise registration of N input images is then posed as the problem of registering N vector images.

An external file that holds a picture, illustration, etc.
Object name is nihms-172754-f0001.jpg

Fig. 1

An attribute vector is the signature of a given voxel (as ^{). The local pattern of an attribute is estimated from the neighborhood centered at location}x, and described in the probabilistic fashion (as ^).

It is worth noting that different attributes should ideally be orthogonal or statistically independent. For example, if the intensity gradients are calculated in median difference manner, the obtained gradients would be dependent only on the two adjacent voxels but independent of the center location under consideration. Then the correlation between any two different attributes approaches to zero. This requirement or assumption simplifies the optimization in registration, in that it enables us to evaluate the contribution of each attribute independently. Further, different image attributes should have different significances according to the specific ongoing stage of registration. In the next section we will explain how to specify the importance of a certain attribute.

2.2 Local Pattern Matching

For intensity guided registration, the common criterion to determine whether registration should be further continued or have already converged is based on the estimation of intensity similarity/dissimilarity, though the metric itself can be formulated in a variety of ways. In the case of groupwise registration, the variance of intensities can be computed for each location, where voxels from different images form a stack. The integration of the variances across the whole image domain results in a cost function, which can be minimized for groupwise registration. However, such a measurement is especially susceptible to imaging distortions and anatomical ambiguities, since intensity from a single voxel is hardly robust.

An improvement to intensity variance is the stack entropy (Balci et al., 2007; Zöllei et al., 2005) which plots the histogram of intensities in each voxelwise stack and calculates the Shannon entropy for the histogram. Nevertheless, to accurately estimate the entropy is challenging when sample size is not large enough (Paninski, 2003). Meanwhile, sampling strategy is widely applied in medical image registration. In Balci et al. (Balci et al., 2007), a subset of voxels, which are randomly sampled from the domain of input images, is used to steer the estimations of deformation fields. The sampling rate can be very low (e.g., 1% of the volume size or even less). Though this strategy can greatly alleviate computation complexity, it raises the concern related to the accuracy and consistency of registration. In particular, a uniformly distributed voxel subset can hardly well reflect the original image, especially when all sampled voxels are equally treated regardless of their intrinsic significances.

For better registration, we introduce the local pattern matching mechanism into the proposed method. Given the center voxel location in a certain image, its neighborhood is regarded as the local pattern of the center voxel. Registration is then formulated as the problem of finding the transformations which could minimize the inhomogeneity among corresponding local patterns with the same center locations in different images. For each location, a PDF could be estimated to reflect the distribution of a given attribute within the neighborhood (as ⁱⁿFig. 1). Following Parzen windowing in probability estimation, the probabilistic distribution of attribute a^j, centered at location x, can be written as:

\Pr (a^{j} ∣ i, x) = \frac{1}{K} \underset{Δ x}{Σ} G_{σ} (a^{j} - a_{i, x + Δ x}^{j})

(1)

where G_σ is a Gaussian kernel with variance σ, K is the overall number of voxels in the neighborhood, and Δx is the offset from the center location. The allowable offset in our method is set within a 3×3×3 cubic neighborhood currently.

The PDF here works as the signature of the local pattern for each voxel. Assuming that all input images are already registered, the series of PDFs {Pr(a^j|i, x)|i = 1, ⋯, N} in the stack should be very similar. Therefore, a well-behaved groupwise registration algorithm should give the minimal divergence for the PDF stack which consists of voxels from the same locations of different images. To measure this divergence, a popular information-theoretic metric, called the Jensen-Shannon (JS) divergence (Lin, 1991), can then be employed:

JS (x, j) = ℋ_{p} (Σ_{i = 1}^{N} π_{i} \cdot \Pr (a^{j} ∣ i, x)) - Σ_{i = 1}^{N} π_{i} \cdot ℋ_{p} (\Pr (a^{j} ∣ i, x))

(2)

where H_p(·) indicates the Shannon entropy, and each individual image is given a factor π_i to weight its contribution. It is worth noting that the JS divergence was adopted in Wang et al. (Wang et al., 2008) for groupwise registration of different point sets, by regarding the whole point set as a global pattern represented by Gaussian mixtures. In contrast, we consider local patterns (PDFs of each attribute in the neighborhood) as representations of regional anatomical structures.

By setting the weighting factor π_i of each image to 1/N, the divergence of the local patterns for the neighborhoods centered at x across different images can be rewritten from Eq. 2 as:

JS (x, j) = \frac{1}{N} Σ_{i = 1}^{N} \int \Pr (a^{j} ∣ i, x) \cdot \log \frac{\Pr (a^{j} ∣ i, x)}{Σ_{s = 1}^{N} \Pr (a^{j} ∣ s, x) ∕ N} {da}^{j}

(3)

Further, given the center location of the neighborhood and the specific attribute under consideration, we could define a quantity V(x,j):

V (x, j) = \log \frac{\Pr (a^{j} ∣ i, x)}{Σ_{s = 1}^{N} \Pr (a^{j} ∣ s, x) ∕ N}

(4)

and thus the local JS divergence in Eq. 3 equals to the conditional expectation of V(x,j):

JS(x, j) = E(V(x, j)∣(a^∣i, x))

(5)

The expectation would approach to zero if all patterns had the same PDF.

2.3 Cost Function and Implementation Issues

Across the whole image domain should be the contributions from all locations and different attributes summed up, in order to form a unified metric which is used as a goodness measurement of the groupwise registration. As mentioned above, the correlation between any two selected attributes is assumed to be close to zero. Therefore, contributions from different attributes are calculated individually and then combined linearly. In this study, we formulate the attribute vector by incorporating intensity as well as the magnitude of the intensity gradient. The reason for utilizing the gradient magnitude is that it can provide accurate guidelines to match local boundaries. In the appendix, we use a toy example to demonstrate the power of the attribute vector and to illustrate the motivation for local pattern matching.

It is, however, not clear how each attribute should be weighted according to its significance and to the stages of registration. We could use the previously reported learning method (Wu et al., 2007) to extract a series of attributes, and determine the optimal combination of their weightings through a training procedure prior to registration. For the sake of simplicity, we design a linearly adaptive assignment strategy to weight these two types of attributes and balance between intensity as well as the gradient magnitude. Empirically, intensity based matching is relatively less sensitive to mismatches of boundaries than gradient. Therefore, we set higher intensity weight in the initial stages of registration. As the registration progresses and images become better aligned, the weight given to the gradient magnitude can be increased, in order to refine the alignment of the structural boundaries. Assuming the weighting factor of the j-th attribute is W^j (constrained by Σ_jW^j = 1), then the overall cost function for minimization can be written as:

\begin{matrix} f & = E (JS (x, j) ∣ x, j) \\ = \int W^{j} \cdot \Pr (x) Σ_{j = 1}^{D} JS (x, j) dx \end{matrix}

(6)

In Eq. 6, the term Pr(x) serves as spatial weighting factors for different locations in the image domain. Traditionally, a subset of voxels is uniformly sampled from the image domain, and a specific weighting factor can be assigned to each voxel location. For example, the optimal scale saliency measurement (Kadir and Brady, 2001) was introduced to mutual information based registration and succeeded in improving the robustness and consistency of the algorithm (Luan et al., 2008). However, the relationship between the cost function and the weighting method is still under investigation. And usually the calculation of potential weighting factors (e.g., optimal scale saliency) is very computation intensive.

To avoid the above problems, we have applied the concept of importance sampling strategy (Bhagalia et al., 2009; Sundar et al., 2007) to the proposed groupwise registration method. Taking advantage of the intensities and the gradient magnitudes already computed for the purpose of registration, we could average, smooth, and then normalize the gradient magnitude values on the overall image domain. The produced values signify the importance of each voxel. Based on the generated importance map, a number of voxels can be sampled to drive the registration. Higher importance implies that the location under consideration is more likely to be drawn in non-uniform sampling, while lower value denotes higher possibility of rejection. Based on the importance map in Fig. 2(b), for instance, an initial subset of voxels, shown in red in Fig. 2(a), can be sampled to guide the registration. Note, with the progress of registration, more and more sampled voxels will gradually participate (red followed by green, and blue) in the registration. The spatial distribution of sampled voxels confirms that samples are more concentrated in boundary abundant areas, where importance values are generally higher.

An external file that holds a picture, illustration, etc.
Object name is nihms-172754-f0002.jpg

Fig. 2

Non-uniform sampling is applied, and a subset of voxel locations is drawn and overlaid on the original slice in (a). Initially, only the red voxels in (a) will be sampled to drive the registration. As registration progresses, more and more voxels indicated in green and blue, will gradually take part in the registration. It is worth noting that most sampled voxels are located at the boundaries of prominent anatomical structures, as indicated by the importance map for the non-uniform sampling, as shown in (b).

Given a subset of voxel locations {x_m, m = 1, ⋯, M} drawn from the non-uniform sampling procedure, the term Pr(x) in Eq. 6 can then be implicitly represented by the intrinsic spatial distribution within the subset. As a result, the overall cost function in Eq. 6 could be expressed as:

f = W^{j} \cdot Σ_{m = 1}^{M} Σ_{j = 1}^{D} JS (x_{m}, j)

(7)

where {x_m, m = 1, ⋯, M} complies with the distribution Pr(x) in Eq. 6.

To minimize the overall f in Eq. 7, we can analytically solve for its gradient and feed it as the steepest descent direction to a first order optimizer. As mentioned above, W^j is associated with the stage of registration. Given a certain iteration in optimization, the configuration of W^j is constant. Therefore, the partial derivative of the cost function can be expressed as:

\frac{\partial f}{\partial T_{i}} = W^{j} \cdot Σ_{m = 1}^{M} Σ_{j = 1}^{D} \frac{\partial JS (x_{m}, j)}{\partial T_{i}}

(8)

where T_i indicates the transformation estimated for the i-th image.

To further improve the robustness, a multi-resolution registration strategy is adopted. In each resolution, except for the highest one, images are down-sampled firstly. The attributes are then calculated on those down-sampled images. This approach ensures that the generated attribute vector is scale-related, thus can better capture anatomical structures in different resolutions. B-Spline transformation model is employed to describe the deformation for each image. To better avoid local minima in registration, the number of B-Spline control points is associated to the resolution of registration. At a coarser resolution, the optimization problem is effectively reduced to a lower degree-of-freedom (DOF) problem, which is achieved by using a smaller number of B-Spline control points. The estimated transformations will be progressively refined as the registration progresses by gradually increasing the number of B-Spline control points.

2.1 Attribute and Attribute Vector

Fig. 1

2.2 Local Pattern Matching

\Pr (a^{j} ∣ i, x) = \frac{1}{K} \underset{Δ x}{Σ} G_{σ} (a^{j} - a_{i, x + Δ x}^{j})

(1)

JS (x, j) = ℋ_{p} (Σ_{i = 1}^{N} π_{i} \cdot \Pr (a^{j} ∣ i, x)) - Σ_{i = 1}^{N} π_{i} \cdot ℋ_{p} (\Pr (a^{j} ∣ i, x))

(2)

By setting the weighting factor π_i of each image to 1/N, the divergence of the local patterns for the neighborhoods centered at x across different images can be rewritten from Eq. 2 as:

JS (x, j) = \frac{1}{N} Σ_{i = 1}^{N} \int \Pr (a^{j} ∣ i, x) \cdot \log \frac{\Pr (a^{j} ∣ i, x)}{Σ_{s = 1}^{N} \Pr (a^{j} ∣ s, x) ∕ N} {da}^{j}

(3)

Further, given the center location of the neighborhood and the specific attribute under consideration, we could define a quantity V(x,j):

V (x, j) = \log \frac{\Pr (a^{j} ∣ i, x)}{Σ_{s = 1}^{N} \Pr (a^{j} ∣ s, x) ∕ N}

(4)

and thus the local JS divergence in Eq. 3 equals to the conditional expectation of V(x,j):

JS(x, j) = E(V(x, j)∣(a^∣i, x))

(5)

The expectation would approach to zero if all patterns had the same PDF.

2.3 Cost Function and Implementation Issues

\begin{matrix} f & = E (JS (x, j) ∣ x, j) \\ = \int W^{j} \cdot \Pr (x) Σ_{j = 1}^{D} JS (x, j) dx \end{matrix}

(6)

Fig. 2

f = W^{j} \cdot Σ_{m = 1}^{M} Σ_{j = 1}^{D} JS (x_{m}, j)

(7)

where {x_m, m = 1, ⋯, M} complies with the distribution Pr(x) in Eq. 6.

\frac{\partial f}{\partial T_{i}} = W^{j} \cdot Σ_{m = 1}^{M} Σ_{j = 1}^{D} \frac{\partial JS (x_{m}, j)}{\partial T_{i}}

(8)

where T_i indicates the transformation estimated for the i-th image.

3 Experimental Results

We will demonstrate in this section the performance of the proposed attributed vector guided groupwise registration using experiments on both real and simulated data. All results are compared with the implementation of the congealing method reported in Balci et al. (Balci et al., 2007), which is designated as the “intensity” guided method in the following. For fair comparison, the two methods are constructed using the same framework, with the same optimizer, transformation model, and multi-resolution registration strategy. For real data, the alignment of the NIREP NA0 data by our method was more accurate than the intensity guided method in terms of both intensity residual errors and ROI overlap ratios. Another experiment was performed on simulated atrophic data. And our method again yielded better results in detecting abnormal atrophies in brain images, indicating its potential clinical applications.

3.1 Anatomical Correspondence

We have employed the NIREP NA0 dataset (Christensen et al., 2006) to demonstrate the capability of the attribute vector guided registration algorithm. In the dataset, there are 16 individual brain MR images, each of which has 32 manual anatomical labels. A typical slice from one subject in the dataset is shown in Fig. 3(a), with its manual labels displayed in Fig. 3(b). For the sake of fairness, the same registration settings are applied to the intensity guided and the attribute vector guided methods. For example, both methods use the same number of levels for the multi-resolution registration framework, and the same number of B-Spline control points in each level.

An external file that holds a picture, illustration, etc.
Object name is nihms-172754-f0003.jpg

Fig. 3

A typical slice (a) from the NA0 dataset, and its manual anatomical labels (b): Different colors indicate different manually delineated anatomical structures.

To quantify the performance, we calculate intensity residual errors between each of the 16 aligned images and their group mean. Then, a standard deviation volume of the intensity residual errors for the 16 images is computed. The histograms of the standard deviation volumes yielded by both two methods are shown in Fig. 4(a). The result of our method (blue curve) shows a distribution much more concentrated to the low error range. This implies that more accurate alignment of the 16 images is achieved by the attribute vector guided method compared to the intensity guided method (red curve). The 3-D rendering views of the standard deviation volumes are also provided in Fig. 4(a). It is worth noting that the two views are adjusted to the same contrast. The relatively darker volume associated with the blue curve confirms that our method leads to fewer mismatches for both cortical and subcortical regions. Besides, in Fig. 4(b), we provide typical slices extracted from each of the mean images yielded by both methods, as well as the 3-D renderings of the two images. Similar conclusion can be reached that our method is superior to the intensity guided registration method. Specifically, more anatomical details are preserved and the produced mean image is sharper.

An external file that holds a picture, illustration, etc.
Object name is nihms-172754-f0004.jpg

Fig. 4

Comparison of histograms of standard deviation volumes in (a) and the produced mean images in (b). The 3D rendering views of standard deviation volumes in (a) confirm that the attribute vector guided groupwise registration leads to lower intensity residual errors than the registration guided by intensity alone. In (b), a typical cortical/sub-cortical slice as well as the 3D rending of the mean image yielded by our method (in blue box) show more abundant details than that of the intensity guided method (in red box).

We further evaluate the neuroanatomical consistency in registration with the help of the manually delineated labels. We warp the labels according to the estimated deformations for comparison. Taking the left and the right precentral gyri (L/R-PCGs) as example, a voxel-wise probability map can be estimated by counting for each voxel the total number of registered images which have the specific L/R-PCG labels. If the probability value is 1.0, all 16 subjects after registration have the same L/R-PCG labels at the voxel under consideration. In Fig. 5, we provide 3D renderings of the two probability maps related to the intensity guided method and the attribute vector guided method, respectively. It can be observed that our method on average produces higher overlap probabilities for L/R-PCGs. Areas with significant improvement are highlighted by the green elliptic contours in Fig. 5.

An external file that holds a picture, illustration, etc.
Object name is nihms-172754-f0005.jpg

Fig. 5

L/R-PCG overlap probability maps of the intensity guided method (left) and the attribute vector guided method (right). It can be observed that the latter yields higher label consistency.

Registration consistency is further evaluated on all anatomical labels by quantitatively comparing the overlap ratios. We boost labels in common space by voting out the majority label across all registered images voxel by voxel. After that, we calculate the intersection as well as the union between the voted label and each warped label. The measure utilized here is the Jaccard coefficient, or overlap ratio, which is defined as the division between sizes of the intersection region and the union region. In Fig. 6, we provide the overall overlap ratios, as well as the standard deviations, of all 32 labels averaged from 16 registered images. The results for the intensity guided and the attribute vector guided methods are indicated by red and blue, respectively. Our method shows consistent improvement in all 32 labels in terms of overlap ratio, with an average increase of 5.06%. Additionally, we define the voted label as the positive while its complement the negative. And then, given a warped label as an observation, the specificity and the sensitivity for that observation can be calculated following:

\begin{matrix} Specificity = \frac{| True Negtive |}{| True Negtive | + | False Positive |} \\ Sensitivity = \frac{| True Positive |}{| True Positive | + | False Negtive |} \end{matrix}

(9)

In Fig. 7, the mean specificity and sensitivity for each individual label are provided along with their standard deviations. For all labels, the attribute vector guided method shows improvement of mean values, which implies more reliable registration results. Furthermore, Fig. 8 shows seven subjects (at the same slice) selected from 16 registered images. As indicated by the arrows, images registered by our method (Fig. 8(b)) are much more consistent than those yielded by the intensity guided method (Fig. 8(a)).

An external file that holds a picture, illustration, etc.
Object name is nihms-172754-f0006.jpg

Fig. 6

The overlap ratios of 32 manually labeled ROIs by both the intensity guided and the attribute vector guided groupwise registration methods are shown in red and blue, respectively. Besides the mean ratio, the standard deviation is also plotted as error-bar for each label. Our method produces better label overlaps for all 32 labels and achieves more consistent registration results.

An external file that holds a picture, illustration, etc.
Object name is nihms-172754-f0007.jpg

Fig. 7

For each label, the specificity and the sensitivity are plotted in left panel and right panel, respectively. For each label, values for both the left (top) and right (bottom) hemispheres are plotted. The cross symbol indicates mean specificity/sensitivity, and the length of the bar shows the standard deviation.

An external file that holds a picture, illustration, etc.
Object name is nihms-172754-f0008.jpg

Fig. 8

Seven selected subjects are shown at the same slice after registration via (a) the intensity guided method and (b) the attribute vector guided method, respectively. Ventricles in row (b) are more similar to each other than those in row (a), especially for the locations indicated by arrows.

3.2 Application: Atrophy Detection

After groupwise registration, all input images are expected to be warped to a common space and their group mean can be taken as an atlas. The morphological diversity in the intensity space is then conveyed by the collection of transformation fields which are estimated for individual subjects. The deformation fields can be used not only to define the anatomical correspondence but also to capture morphological differences between input image and the atlas (Baloch et al., 2007; Davatzikos et al., 1996). For example, using the estimated deformation fields, quantitative analysis can be performed to characterize the pathological abnormalities between normal controls and dementia subjects (Christensen et al., 1993; John and Karl, 1999; Miller et al., 1997; Miller and Younes, 2001; Shen and Davatzikos, 2002).

In this experiment, we use two sets of images, each of which contains 12 images. The first set is composed of 12 images as control subjects, and the second set is obtained from the first set by introducing simulated atrophies in the neighborhoods of both the pre-central gyrus (PCG) and the superior temporal gyrus (STG) (Davatzikos et al., 2001; Xue et al., 2006). Sample slices are shown in Fig. 9(a) and (b), where atrophies in PCG and STG are highlighted. Also, the difference between normal control and atrophic subject is shown at the bottom of Fig. 9. Note that the contrast of the differences has been adjusted for better visual inspection. The amount of atrophy, in both STG and PCG respectively, is around 10% reflecting realistic abnormalities in aging or pathological brains (Desikan et al., 2009; Ferdinando et al., 2009; Petersen et al., 1999; Petersen et al., 2005).

An external file that holds a picture, illustration, etc.
Object name is nihms-172754-f0009.jpg

Fig. 9

Control images, and simulated images with atrophy introduced to the pre-central gyrus (PCG) in (a) and the superior temporal gyrus (STG) in (b). The difference between the control image and the atrophic image is provided in the bottom. Note that the contrast of the two differences has been adjusted, for better visual inspection.

After groupwise registration, the Jacobian determinant map corresponding to each deformation is computed. The Jacobian determinant value measures the change of the voxel volume. Considering the fact that the determinant equal to 1.0 represents a voxelwise identical transformation, a >1.0 determinant represents expansion of the volume while a <1.0 determinant implies shrinkage. Based on the Jacobian determinant maps, we perform analysis on scalar values, in contrast to the dense deformation fields which are described by 3D vectors.

With the software package SPM (www.fil.ion.ucl.ac.uk/spm/), a two-sample t-test is performed on the Jacobian determinant maps to compare the group differences between the two sets of images. The two sample t-test can identify significant group differences and ignore the random distortions. Under the same confidence interval, a higher t-value indicates a more salient group difference occurring at that specific location. Detailed results of the t-test are provided in Table 1, which confirms that our method yields improved t-values by 12.95% in PCG and 13.99% in STG. It implies that the estimated deformations of our registration method can more accurately capture subtle morphological changes (e.g. tiny atrophies) within input images. This improvement is crucial for automatic atrophy detection in potential clinical applications.

Table 1

Two sample t-test results of simulated brain atrophy detection.

		Intensity Guided			Attribute Vector Guided
		p_FWE	p_FDR	t-value	p_FWE	p_FDR	t-value

Locations	PCG	0.043	0.006	*6.72*	0.016	0.003	*7.59*

	STG	0.016	0.006	*7.29*	0.005	0.003	*8.31*

3.1 Anatomical Correspondence

Fig. 3

A typical slice (a) from the NA0 dataset, and its manual anatomical labels (b): Different colors indicate different manually delineated anatomical structures.

Fig. 4

Fig. 5

L/R-PCG overlap probability maps of the intensity guided method (left) and the attribute vector guided method (right). It can be observed that the latter yields higher label consistency.

\begin{matrix} Specificity = \frac{| True Negtive |}{| True Negtive | + | False Positive |} \\ Sensitivity = \frac{| True Positive |}{| True Positive | + | False Negtive |} \end{matrix}

(9)

Fig. 6

Fig. 7

Fig. 8

3.2 Application: Atrophy Detection

Fig. 9

Table 1

Two sample t-test results of simulated brain atrophy detection.

		Intensity Guided			Attribute Vector Guided
		p_FWE	p_FDR	t-value	p_FWE	p_FDR	t-value

Locations	PCG	0.043	0.006	*6.72*	0.016	0.003	*7.59*

	STG	0.016	0.006	*7.29*	0.005	0.003	*8.31*

4. Conclusion and Discussion

We have developed a novel attribute vector guided groupwise registration algorithm which is guided by attribute vectors. Although intensity similarity is a key factor and important in image registration, integrating more attributes into registration has been proven useful. In our algorithm, different attributes which capture different aspects of regional characteristics are used as structural signature for each voxel. Registration is then achieved by using a local pattern matching mechanism and by minimizing the attribute PDF dissimilarity across all subjects. Contributions from different attributes and different locations are adaptively weighted, and are mathematically integrated into a cost function which is minimized by following its steepest descent direction. More consistent and accurate registration results can thus be achieved, as experimental results confirm less residual intensity error and higher overlap ratio. Our future work includes incorporating more image attributes to further improve the groupwise registration accuracy. Also, we will apply our method to real brain data to seek for its potential clinical usages.

^{Department of Computer Science University of North Carolina at Chapel Hill, NC 27599}ude.cnu.sc@gnawnaiq

^{Department of Radiology and BRIC University of North Carolina at Chapel Hill, NC 27599}ude.cnu.dem@uwrg, ude.cnu.dem@paytp, ude.cnu.dem@nehsgd

Publisher's Disclaimer

Abstract

Groupwise registration has been recently introduced to simultaneously register a group of images by avoiding the selection of a particular template. To achieve this, several methods have been proposed to take advantage of information-theoretic entropy measures based on image intensity. However, simplistic utilization of voxelwise image intensity is not sufficient to establish reliable correspondences, since it lacks important contextual information. Therefore, we explore the notion of attribute vector as the voxel signature, instead of image intensity, to guide the correspondence detection in groupwise registration. In particular, for each voxel, the attribute vector is computed from its multi-scale neighborhoods, in order to capture the geometric information at different scales. The probability density function (PDF) of each element in the attribute vector is then estimated from the local neighborhood, providing a statistical summary of the underlying anatomical structure in that local pattern. Eventually, with the help of Jensen-Shannon (JS) divergence, a group of subjects can be aligned simultaneously by minimizing the sum of JS divergences across the image domain and all attributes. We have employed our groupwise registration algorithm on both real (NIREP NA0 dataset) and simulated data (12 pairs of normal control and simulated atrophic dataset). The experimental results demonstrate that our method yields better registration accuracy, compared with a popular groupwise registration method.

Keywords: Groupwise registration, attribute vector, Jensen-Shannon divergence

Abstract

Footnotes

Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errorsmaybe discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Footnotes

Attribute Vector Guided Groupwise Registration

1. Introduction

2. Methods

2.1 Attribute and Attribute Vector

2.2 Local Pattern Matching

2.3 Cost Function and Implementation Issues

2.1 Attribute and Attribute Vector

2.2 Local Pattern Matching

2.3 Cost Function and Implementation Issues

3 Experimental Results

3.1 Anatomical Correspondence

3.2 Application: Atrophy Detection

Table 1

3.1 Anatomical Correspondence

3.2 Application: Atrophy Detection

Table 1

4. Conclusion and Discussion

Abstract

Footnotes

Reference