Difference between revisions of "Projects:GroupwiseRegistration"

From NAMIC Wiki
Jump to: navigation, search
 
(53 intermediate revisions by 7 users not shown)
Line 1: Line 1:
Back to [[NA-MIC_Collaborations|NA-MIC_Collaborations]], [[Algorithm:MIT|MIT Algorithms]]
+
Back to [[NA-MIC_Internal_Collaborations:StructuralImageAnalysis|NA-MIC Collaborations]], [[Algorithm:MIT|MIT Algorithms]]
 
+
__NOTOC__
 
= Non-rigid Groupwise Registration =
 
= Non-rigid Groupwise Registration =
  
In this work,
+
We aim at providing efficient groupwise registration algorithms
we extend a previously demonstrated entropy based groupwise registration method  
+
for population analysis of anatomical structures.
 +
Here we extend a previously demonstrated entropy based groupwise registration method  
 
to include a free-form deformation model based on B-splines.  
 
to include a free-form deformation model based on B-splines.  
 
We provide  
 
We provide  
Line 17: Line 18:
 
= Description =
 
= Description =
  
[[Image:GroupwiseBspline.png|thumb|200px|Figure 1: desc]]
+
We first describe
[[Image:GroupwiseStackBiModal.PNG|thumb|200px|Figure 1: desc]]
+
the stack entropy cost function and the B-spline based deformation model.
[[Image:GroupwiseMultiResolution.PNG|thumb|200px|Figure 1: desc]]
+
Then we discuss implementation details.  
[[Image:GroupwiseIncreasingScale.PNG|thumb|200px|Figure 1: desc]]
+
Next, we compare groupwise registration to the pairwise method and
[[Image:GroupwiseMeanImages.png|thumb|200px|Figure 1: desc]]
+
evaluate both methods using label prediction values.
[[Image:GroupwiseBarsWMGM.png|thumb|200px|Figure 1: desc]]
+
 
[[Image:GroupwiseBarsManual.png|thumb|200px|Figure 1: desc]]
 
  
'''Objective Function'''
+
''Objective Function''
  
 +
[[Image:GroupwiseStackBiModal.PNG|thumb|350px|Figure 1: On the left is shown a stack of images
 +
and a sample pixel stack around a cortical region. On the left is shown the Gaussian(blue) fittet to
 +
a real sample from the dataset we used along with the non-parametric density estimate(red).
 +
Note that the distribution is bi-modal because of white matter-gray matter transaction.]]
  
 
In order to align all subjects in the population,
 
In order to align all subjects in the population,
Line 36: Line 40:
 
subjects are simultenously driven to the common tendency of the population.
 
subjects are simultenously driven to the common tendency of the population.
  
We employ a Parzen window based density estimation scheme to estimate univariate entropies \cite{duda}:
+
We employ a kernel based density estimation scheme to estimate univariate entropies.
\begin{equation}
 
f = -\sum_{v=1}^{V} \frac{1}{N} \sum_{i=1}^{N} \log \frac{1}{N}\sum_{j=1}^{N} G_\sigma(d_{ij}(x_v) )
 
\end{equation}
 
where $d_{ij}(x)=I_i(T_i(x))-I_j(T_j(x))$ is the distance between intensity values
 
of a pair of images evaluated at a point in the reference frame
 
and $G_\sigma$ is a Gaussian kernel with variance $\sigma^2$.
 
The objective function achieves its minimum when the intensity differences are small.  
 
 
Using the entropy measure we obtain a better treatment of transitions between different
 
Using the entropy measure we obtain a better treatment of transitions between different
 
tissue types, such as gray matter-white matter transitions in the cortical regions
 
tissue types, such as gray matter-white matter transitions in the cortical regions
where intensity distributions can be multi-modal.
+
where intensity distributions can be bi-modal as shown in Figure 1.
Parzen window density estimator allows us to obtain
+
 
analytical expressions for the derivative of the objective function
 
with respect to transform parameters
 
  
'''Deformation Model'''
+
 
 +
''Deformation Model''
 +
 
 +
[[Image:GroupwiseBspline.png|thumb|350px|Figure 2: An example deformation field. The local neighborhood affecting the deformation is overlayed on the image.]]
  
 
For the nonrigid deformation model,
 
For the nonrigid deformation model,
 
we define a combined transformation consisting of  
 
we define a combined transformation consisting of  
 
a global and a local component
 
a global and a local component
\begin{equation}
+
 
 +
:<math>
 
T(\mathbf{x}) = T_{local}({T_{global}(\mathbf{x})})
 
T(\mathbf{x}) = T_{local}({T_{global}(\mathbf{x})})
\end{equation}
+
</math>
where $T_{global}$ is a twelve parameter affine transform and
 
$T_{local}$ is a deformation model based on B-splines.
 
  
Following Rueckert et al.'s formulation \cite{rueckert},
+
where <math>T_{global}</math> is a twelve parameter affine transform and
we let $\mathbf{\Phi}$ denote an
+
<math>T_{local}</math> is a deformation model based on B-splines.
$n_x \times n_y \times n_z $ grid of control points $\Phi_{i,j,k}$ with uniform
+
 
spacing. The free form deformation can be written as the 3-D tensor product
+
The free form deformation can be written as the 3-D tensor product
 
of 1-D cubic B-splines.
 
of 1-D cubic B-splines.
\begin{equation}
+
 
 +
:<math>
 
T_{local}(\mathbf{x}) = \mathbf{x} + \sum_{l=0}^3\sum_{m=0}^3\sum_{n=0}^3 B_l(u)B_m(v)B_n(w) \Phi_{i+l,j+m,k+n}
 
T_{local}(\mathbf{x}) = \mathbf{x} + \sum_{l=0}^3\sum_{m=0}^3\sum_{n=0}^3 B_l(u)B_m(v)B_n(w) \Phi_{i+l,j+m,k+n}
\end{equation}
+
</math>
where
 
$\mathbf{x}=(x,y,z)$,
 
$i=\lfloor x/n_x \rfloor  -1$,
 
$j=\lfloor y/n_y\rfloor -1$,
 
$k=\lfloor z/n_z\rfloor -1$,
 
$u=x/n_x -\lfloor x/n_x\rfloor$,
 
$v=y/n_y -\lfloor y/n_y\rfloor$,
 
$w=z/n_z -\lfloor z/n_z\rfloor$
 
and where $B_l$ is $l$'th cubic B-spline basis function.
 
Using the same expressions for $u,v$ and $w$ as above, the derivative of the deformation field with respect to B-spline coefficients
 
can be given by
 
\begin{equation}
 
\frac{\partial T_{local}(x,y,z) }{\partial \Phi_{i,j,k}} = B_l(u)B_m(v)B_n(w)
 
\end{equation}
 
where $l=i-\lfloor x/n_x\rfloor+1$, $m=j-\lfloor y/n_y\rfloor+1$ and $n=k-\lfloor z/n_z\rfloor+1$. We consider $B_l(u) = 0$ for $l<0 $ and $l>3$.
 
The derivative terms are nonzero only in the neighborhood of a given point. Therefore,
 
optimization of the objective function using gradient descent can be implemented efficiently.
 
  
  
%The resolution of the deformation field
+
%can be controlled by the grid size.
+
where <math>B_l</math> is <math>l</math>'th cubic B-spline basis function. <math>(u,v,w)</math> is the distance
%In order to capture shape variations at different resolution levels
+
to <math>(x,y,z)</math> from the control point <math>\Phi_{i,j,k}</math> as shown in Figure 2.
%we start the registration at a low resolution level with a coarse grid
 
%and increase the resolution by refining the grid of control points.
 
%Results of a registration at a coarse level are used to initialize
 
%the grid of control points at a higher resolution level.
 
  
As none of the images are chosen as an anatomical reference,
+
The deformation of a given point can be found using only the control points in the neighborhood of the given point. Therefore,
it is necessary to add a geometric constraint to define the reference coordinate frame.
+
optimization of the objective function can be implemented efficiently.
Similar to Bhatia et al. \cite{bhatia}, we define the reference frame
 
by constraining the average deformation to be the identity transform:
 
\begin{equation}\label{eq:1}
 
\frac{1}{N}\sum_{n=1}^{N} T_n(\mathbf{x}) = \mathbf{x}
 
\end{equation}
 
  
This constraint assures that the reference frame lies in the center of the
 
population. In the case of B-splines, the constraint can be satisfied
 
by constraining the sum of B-spline coefficients across images
 
to be zero. In the gradient descent optimization scheme,
 
the constraint can be forced by
 
subtracting the mean from each update vector \cite{bhatia}.
 
  
 +
''Implementation''
  
 
+
[[Image:GroupwiseIncreasingScale.PNG|thumb|350px|Figure 3: A registration schedule using gradually increasing deformation field complexity. From left to right deformation fields for increasing deformation field complexity. ]]
We compare our groupwise algorithm to a pairwise method where we register
 
each subject to the mean intensity using sum of square differences.
 
The objective function for pairwise registration to the mean can be described as follows
 
\begin{equation}
 
f_{pair} = \sum_{n=1}^{N} ( I_n(T_n(x)) - \mu(x)  )^2
 
\end{equation}
 
where $\mu$ is defined as the mean of the intensities $\mu(x) = \frac{1}{N}\sum_{n=1}^{N} I_n(T_n(x))$.
 
During each iteration we consider the mean image as a reference image
 
and register
 
every subject to the mean image using sum of squared differences.  
 
After each iteration the mean image is updated and pairwise registrations are performed until convergence.
 
 
 
 
 
 
 
The images in Figure \ref{fig:mean} show central slices of 3D images after registration.
 
Visually, mean images get sharper and variance images becomes darker, especially around central ventricles and cortical regions.
 
We can observe that anatomical variability at cortical regions causes significant blur for
 
GM, WM and CSF structures using affine registration.  
 
Finer scales of B-spline deformation fields capture a significant part of this anatomical variability and
 
the tissue label overlap images get sharper.
 
 
 
 
 
 
 
We evaluate registration results by measuring label prediction accuracy
 
in a leave-one-out fashion for the two different sets of segmentation labels.
 
To predict the segmentation labels of a subject, we
 
use majority voting for the labels in the rest of the population.
 
We compute the Dice measure between the predicted and the true labels and average
 
over the whole population.
 
%Table \ref{table:pred} shows average Dice measures for two automatically segmented
 
%tissue types: grey matter (GM) and white matter (WM) in the first two columns.
 
%We excluded cerebro-spinal fluid (CSF) from the evaluation because CSF labels
 
%near the cortex get corrupted
 
%after preprocessing the data with skull stripping.
 
%In the last two columns of Table \ref{table:pred} we display average prediction values for
 
%manual segmentations of the cortical structures: superior temporal gyrus and para-hippocampus in the
 
%left and right hemispheres; the subcortical structures: amygdala and hippocampus in the left and
 
%right hemispheres.
 
 
 
 
 
 
 
'''Implementation'''
 
  
 
We provide an efficient optimization scheme by using line search with the gradient descent algorithm.
 
We provide an efficient optimization scheme by using line search with the gradient descent algorithm.
For computational efficiency, we employ a stochastic subsampling procedure \cite{pluim}.  
+
For computational efficiency, we employ a stochastic subsampling procedure.  
 
In each iteration of the algorithm,  
 
In each iteration of the algorithm,  
 
a random subset is drawn from all samples and the objective function is evaluated
 
a random subset is drawn from all samples and the objective function is evaluated
 
only on this sample set.  
 
only on this sample set.  
The number of samples to be used in this method depends on the number of the parameters of
 
the deformation field.
 
Using the number of samples on the order of the number of variables works well in practice.
 
  
 
To obtain a dense deformation field capturing anatomical variations at different scales,
 
To obtain a dense deformation field capturing anatomical variations at different scales,
 
we gradually increase the complexity of the deformation field by refining the grid of B-spline control points.
 
we gradually increase the complexity of the deformation field by refining the grid of B-spline control points.
First, we perform a global registration with affine transforms.
+
 
Then, we use affine transforms to initialize a low resolution deformation field
+
 
at a grid spacing around $32$~voxels.
+
[[Image:GroupwiseMultiResolution.PNG|thumb|350px|Figure 4: An example showing the multi-resolution scheme. The registration is first performed at a coarse scale by downsampling the input.
We increase the resolution of the deformation field to $8$~voxels
+
Results from coarser scales are used to initialize
by using registration results at coarser grids to initialize finer grids.
+
optimization at finer scales. Also note that the objective function is only evaluated on a small subset of input points. ]]
  
 
As in every iterative search algorithm, local minima pose a significant problem.  
 
As in every iterative search algorithm, local minima pose a significant problem.  
 
To avoid local minima we use a multi-resolution optimization scheme for each resolution level of the deformation field.
 
To avoid local minima we use a multi-resolution optimization scheme for each resolution level of the deformation field.
The registration is first performed at a coarse scale by downsampling the input.
 
Results from coarser scales are used to initialize
 
optimization at finer scales.
 
For each resolution level of the deformation field we used a multi-resolution scheme of three
 
image resolution levels.
 
  
We implemented our groupwise registration method in a multi-threaded fashion using Insight Toolkit (ITK)
+
We implemented our groupwise registration method in a multi-threaded fashion using Insight Toolkit(ITK)
and made the implementation publicly available \cite{itk}.  
+
and made the implementation publicly available [http://www.na-mic.org/svn/NAMICSandBox/trunk/MultiImageRegistration/ (code)].
We run experiments using a dataset of 50 MR images
 
with $256\times256\times128$ voxels on a workstation with four CPUs and
 
8GB of memory. The running time of the algorithm
 
is about 20 minutes for affine registration and two days for non-rigid registration of the entire dataset.
 
The memory requirement of the algorithm depends linearly on the number of input images
 
and was around 3GB for our dataset.
 
  
 +
''Results''
  
 +
[[Image:GroupwiseMeanImages.png|thumb|350px|Figure 5:  Central slices of 3D volumes for groupwise registration. Rows show mean and standard deviation images followed by label overlap images for GM, WM and CSF labels. Columns display the results for affine and B-splines with grid spacing 32, 16 and 8 voxels, respectively. ]]
 +
[[Image:GroupwiseBarsWMGM.png|thumb|350px|Figure 6: GM, WM DICE measures computed for different deformation field resolution levels. Blue bars show the results for groupwise registration and the red bars show the results for registration to the mean setting.]]
 +
[[Image:GroupwiseBarsManual.png|thumb|350px|Figure 7: DICE measures for manually segmented labels. Bars correspond to the same setting as in figure 6. ]]
  
'''Progress'''
 
  
We tested the groupwise registration algorithm on a MR brain dataset and
+
We tested the groupwise registration algorithm on a MR brain dataset.  
compared the results to a pairwise approach \cite{joshi}.  
 
 
The dataset consists of 50 MR brain images of three subgroups:
 
The dataset consists of 50 MR brain images of three subgroups:
 
schizophrenics, affected disorder and normal control patients.  
 
schizophrenics, affected disorder and normal control patients.  
MR images are T1 scans with $256\times256\times128$ voxels  
+
MR images are T1 scans with 256x256x128 voxels  
and $0.9375\times0.9375\times1.5$~$\mbox{mm}^3$ spacing.
+
and 0.9375x0.9375x1.5 mm<sup>3</sup> spacing.  
MR images are preprocessed by skull stripping.
+
For each image in the dataset, an automatic tissue classification
To account for global intensity variations, images are normalized
 
to have zero mean and unit variance.
 
For each image in the dataset, an automatic tissue classification \cite{pohl}
 
 
was performed, yielding gray matter (GM), white matter (WM) and cerebro-spinal
 
was performed, yielding gray matter (GM), white matter (WM) and cerebro-spinal
 
fluid (CSF) labels. In addition, manual segmentations of four subcortical regions  
 
fluid (CSF) labels. In addition, manual segmentations of four subcortical regions  
Line 210: Line 120:
 
gyrus and para-hippocampus) were available for each MR image.
 
gyrus and para-hippocampus) were available for each MR image.
  
 
The prediction accuracy reported in Table \ref{table:pred} is lower than what is typically
 
achieved by segmentation methods. This is to be expected as our
 
relatively simple label prediction method only considers voxelwise
 
majority of the labels in the population and does not use the novel
 
image intensity to predict the labels. Effectively, Table \ref{table:pred} reports
 
the accuracy of the spatial prior (atlas) in predicting the labels
 
before the local intensity is used to refine the segmentation.
 
 
Increasing the complexity of the deformation model improves the
 
Increasing the complexity of the deformation model improves the
 
accuracy of prediction. An interesting open problem is automatically
 
accuracy of prediction. An interesting open problem is automatically
Line 237: Line 139:
 
Comparing the groupwise registration to the pairwise approach, we
 
Comparing the groupwise registration to the pairwise approach, we
 
observe that the sharpness of the mean images and the tissue overlaps
 
observe that the sharpness of the mean images and the tissue overlaps
in Figure~\ref{fig:mean} look visually similar. From Table~\ref{table:pred}, we note that
+
in Figure 5 look visually similar. From Figures 6 and 7, we note that
 
groupwise registration performs slightly better than the pairwise
 
groupwise registration performs slightly better than the pairwise
 
setting in most of the cases, especially as we increase the complexity
 
setting in most of the cases, especially as we increase the complexity
Line 246: Line 148:
 
comparative studies are needed of the two approaches.
 
comparative studies are needed of the two approaches.
  
%Average Dice measures in Table \ref{table:pred} indicate relatively low values, especially for
+
We compare our groupwise algorithm to a pairwise method where we register
%manual segmentations of cortical structures.
+
each subject to the mean intensity using sum of square differences.
%This is to be expected as our relatively simple label prediction method
+
During each iteration we consider the mean image as a reference image
%only considers voxelwise majority of the labels in the population.
+
and register
%Higher overlap values can be obtained by using more advanced segmentation techniques which make us of the
+
every subject to the mean image using sum of squared differences.
% shapes of the structures.
+
After each iteration the mean image is updated and pairwise registrations are performed until convergence.
%As our method keeps the voxelwise intensity distributions in the population, these distributions
 
%can be used to supply prior information to a more advanced segmentation algorithm to increase
 
%the segmentation accuracy.
 
  
 +
The images in Figure 5 show central slices of 3D images after registration.
 +
Visually, mean images get sharper and variance images becomes darker, especially around central ventricles and cortical regions.
 +
We can observe that anatomical variability at cortical regions causes significant blur for
 +
GM, WM and CSF structures using affine registration.
 +
Finer scales of B-spline deformation fields capture a significant part of this anatomical variability and
 +
the tissue label overlap images get sharper.
  
%Evaluation of a groupwise registration
+
=Asymmetric Image-Template Registration=
%challenging task as it is hard to distinguish between registration inaccuracy and anatomical variability
 
%by considering Dice measure as the only evaluation criterion.
 
%However, comparing overlap measures relative to each other gives us useful information.
 
%The blurriness of the mean images in figure \ref{fig:mean} indicate high anatomical variabilities
 
%at cortical structures.
 
%Part of this anatomical variability is captured by higher resolution levels of B-splines,
 
%as can be observed from the increase in Dice measures for cortical structures in Table \ref{table:pred}.
 
%White matter and grey matter tissue type have a significant component in cortical regions; therefore,
 
%the overlap measures for these tissue types also increase as the resolution of the deformation field increases.
 
%Label overlap values for subcortical structures do not improve significantly with an increase in the resolution of the
 
%deformation field.
 
%These structures do not highly correlate with local intensity variations; therefore,
 
%we believe that
 
%the low Dice measures are mostly due to anatomical variabilities in these structures.
 
  
 +
A natural requirement in pairwise image registration is that the resulting deformation is independent of the order of the images. This constraint is typically achieved via a symmetric cost function and has been shown to reduce the effects of local optima. Consequently, symmetric registration has been successfully applied to pairwise image registration as well as the spatial alignment of individual images with a template. However, recent work has shown that the relationship between
 +
an image and a template is fundamentally asymmetric. In this work, we develop a method that reconciles the practical advantages of symmetric registration with the asymmetric nature of image-template registration by adding a simple correction factor to the symmetric cost function. We instantiate our model within a log-domain diffeomorphic registration
 +
framework. Our experiments show exploiting the asymmetry in image-template registration improves alignment in the image coordinates.
  
 +
= Key Investigators =
  
%As for the comparison of groupwise to pairwise approach,  
+
* MIT: Mert R. Sabuncu, B.T. Thomas Yeo, Koen Van Leemput, Serdar K. Balci, Polina Golland.
%we can observe that the sharpness of the mean images and the tissue overlaps in figure \ref{fig:mean} look visually similar.  
+
* Harvard: Sylvain Bouix, Martha E. Shenton, Bruce Fischl, W.M. (Sandy) Wells.
%From figure \ref{table:pred} we note that groupwise registration performs
+
* Kitware: Brad Davis, Louis Ibanez.
%slightly better than the pairwise setting in most of the cases.  
 
%This suggests that considering the population as a whole and registering subjects jointly brings the population into better %alignment than matching each subject to a mean template image.
 
%A mean template image cannot fully explain multi-modal distributions; our groupwise setting
 
%captures the modes in the distributions by considering stack entropies as an alignment criterion.
 
  
= Key Investigators =
+
= Publications =
  
* MIT: Serdar K Balci, Polina Golland, Sandy Wells
 
  
= Publications =
+
[http://www.na-mic.org/publications/pages/display?search=Projects%3AGroupwiseRegistration&submit=Search&words=all&title=checked&keywords=checked&authors=checked&abstract=checked&sponsors=checked&searchbytag=checked| NA-MIC Publications Database on Groupwise Registration]
  
* S.K. Balci, P. Golland, M.E. Shenton, W.M. Wells III. Free-Form B-spline Deformation Model for Groupwise Registration. In Proceedings of MICCAI 2007 Statistical Registration Workshop: Pair-wise and Group-wise Alignment and Atlas Formation, 23-30, 2007.
+
[[Category: Registration]] [[Category:Segmentation]] [[Category:MRI]] [[Category:Schizophrenia]]

Latest revision as of 20:07, 11 May 2010

Home < Projects:GroupwiseRegistration
Back to NA-MIC Collaborations, MIT Algorithms

Non-rigid Groupwise Registration

We aim at providing efficient groupwise registration algorithms for population analysis of anatomical structures. Here we extend a previously demonstrated entropy based groupwise registration method to include a free-form deformation model based on B-splines. We provide an efficient implementation using stochastic gradient descents in a multi-resolution setting. We demonstrate the method in application to a set of 50 MRI brain scans and compare the results to a pairwise approach using segmentation labels to evaluate the quality of alignment. Our results indicate that increasing the complexity of the deformation model improves registration accuracy significantly, especially at cortical regions.

Description

We first describe the stack entropy cost function and the B-spline based deformation model. Then we discuss implementation details. Next, we compare groupwise registration to the pairwise method and evaluate both methods using label prediction values.


Objective Function

Figure 1: On the left is shown a stack of images and a sample pixel stack around a cortical region. On the left is shown the Gaussian(blue) fittet to a real sample from the dataset we used along with the non-parametric density estimate(red). Note that the distribution is bi-modal because of white matter-gray matter transaction.

In order to align all subjects in the population, we consider sum of pixelwise entropies as a joint alignment criterion. The justification for this approach is that if the images are aligned properly, intensity values at corresponding coordinate locations from all the images will form a low entropy distribution. This approach does not require the use of a reference subject; all subjects are simultenously driven to the common tendency of the population.

We employ a kernel based density estimation scheme to estimate univariate entropies. Using the entropy measure we obtain a better treatment of transitions between different tissue types, such as gray matter-white matter transitions in the cortical regions where intensity distributions can be bi-modal as shown in Figure 1.


Deformation Model

Figure 2: An example deformation field. The local neighborhood affecting the deformation is overlayed on the image.

For the nonrigid deformation model, we define a combined transformation consisting of a global and a local component

[math] T(\mathbf{x}) = T_{local}({T_{global}(\mathbf{x})}) [/math]

where [math]T_{global}[/math] is a twelve parameter affine transform and [math]T_{local}[/math] is a deformation model based on B-splines.

The free form deformation can be written as the 3-D tensor product of 1-D cubic B-splines.

[math] T_{local}(\mathbf{x}) = \mathbf{x} + \sum_{l=0}^3\sum_{m=0}^3\sum_{n=0}^3 B_l(u)B_m(v)B_n(w) \Phi_{i+l,j+m,k+n} [/math]


where [math]B_l[/math] is [math]l[/math]'th cubic B-spline basis function. [math](u,v,w)[/math] is the distance to [math](x,y,z)[/math] from the control point [math]\Phi_{i,j,k}[/math] as shown in Figure 2.

The deformation of a given point can be found using only the control points in the neighborhood of the given point. Therefore, optimization of the objective function can be implemented efficiently.


Implementation

Figure 3: A registration schedule using gradually increasing deformation field complexity. From left to right deformation fields for increasing deformation field complexity.

We provide an efficient optimization scheme by using line search with the gradient descent algorithm. For computational efficiency, we employ a stochastic subsampling procedure. In each iteration of the algorithm, a random subset is drawn from all samples and the objective function is evaluated only on this sample set.

To obtain a dense deformation field capturing anatomical variations at different scales, we gradually increase the complexity of the deformation field by refining the grid of B-spline control points.


Figure 4: An example showing the multi-resolution scheme. The registration is first performed at a coarse scale by downsampling the input. Results from coarser scales are used to initialize optimization at finer scales. Also note that the objective function is only evaluated on a small subset of input points.

As in every iterative search algorithm, local minima pose a significant problem. To avoid local minima we use a multi-resolution optimization scheme for each resolution level of the deformation field.

We implemented our groupwise registration method in a multi-threaded fashion using Insight Toolkit(ITK) and made the implementation publicly available (code).

Results

Figure 5: Central slices of 3D volumes for groupwise registration. Rows show mean and standard deviation images followed by label overlap images for GM, WM and CSF labels. Columns display the results for affine and B-splines with grid spacing 32, 16 and 8 voxels, respectively.
Figure 6: GM, WM DICE measures computed for different deformation field resolution levels. Blue bars show the results for groupwise registration and the red bars show the results for registration to the mean setting.
Figure 7: DICE measures for manually segmented labels. Bars correspond to the same setting as in figure 6.


We tested the groupwise registration algorithm on a MR brain dataset. The dataset consists of 50 MR brain images of three subgroups: schizophrenics, affected disorder and normal control patients. MR images are T1 scans with 256x256x128 voxels and 0.9375x0.9375x1.5 mm3 spacing. For each image in the dataset, an automatic tissue classification was performed, yielding gray matter (GM), white matter (WM) and cerebro-spinal fluid (CSF) labels. In addition, manual segmentations of four subcortical regions (left and right hippocampus and amygdala) and four cortical regions (left and right superior temporal gyrus and para-hippocampus) were available for each MR image.

Increasing the complexity of the deformation model improves the accuracy of prediction. An interesting open problem is automatically identifying the appropriate deformation complexity before the registration overfits and the accuracy of prediction goes down. We also note that the alignment of the subcortical structures is much better than that of the cortical regions. It is not surprising as the registration algorithm does not use the information about geometry of the cortex to optimize the alignment of the cortex. In addition, it has been often observed that the cortical structures exhibit higher variability across subjects when considered in the 3D volume rather than modelled on the surface.

Our experiments highlight the need for further research in developing evaluation criteria for image alignment. We used the standard Dice measure, but it is not clear that this measurement captures all the nuances of the resulting alignment.

Comparing the groupwise registration to the pairwise approach, we observe that the sharpness of the mean images and the tissue overlaps in Figure 5 look visually similar. From Figures 6 and 7, we note that groupwise registration performs slightly better than the pairwise setting in most of the cases, especially as we increase the complexity of the warp. This suggests that considering the population as a whole and registering subjects jointly brings the population into better alignment than matching each subject to a mean template image. However, the advantage shown here is only slight; more comparative studies are needed of the two approaches.

We compare our groupwise algorithm to a pairwise method where we register each subject to the mean intensity using sum of square differences. During each iteration we consider the mean image as a reference image and register every subject to the mean image using sum of squared differences. After each iteration the mean image is updated and pairwise registrations are performed until convergence.

The images in Figure 5 show central slices of 3D images after registration. Visually, mean images get sharper and variance images becomes darker, especially around central ventricles and cortical regions. We can observe that anatomical variability at cortical regions causes significant blur for GM, WM and CSF structures using affine registration. Finer scales of B-spline deformation fields capture a significant part of this anatomical variability and the tissue label overlap images get sharper.

Asymmetric Image-Template Registration

A natural requirement in pairwise image registration is that the resulting deformation is independent of the order of the images. This constraint is typically achieved via a symmetric cost function and has been shown to reduce the effects of local optima. Consequently, symmetric registration has been successfully applied to pairwise image registration as well as the spatial alignment of individual images with a template. However, recent work has shown that the relationship between an image and a template is fundamentally asymmetric. In this work, we develop a method that reconciles the practical advantages of symmetric registration with the asymmetric nature of image-template registration by adding a simple correction factor to the symmetric cost function. We instantiate our model within a log-domain diffeomorphic registration framework. Our experiments show exploiting the asymmetry in image-template registration improves alignment in the image coordinates.

Key Investigators

  • MIT: Mert R. Sabuncu, B.T. Thomas Yeo, Koen Van Leemput, Serdar K. Balci, Polina Golland.
  • Harvard: Sylvain Bouix, Martha E. Shenton, Bruce Fischl, W.M. (Sandy) Wells.
  • Kitware: Brad Davis, Louis Ibanez.

Publications

NA-MIC Publications Database on Groupwise Registration