# Bayesian Spatial Binary Regression for Label Fusion in Structural
Neuroimaging

D. Andrew Brown School of Mathematical and Statistical Sciences, Clemson
University, Clemson, SC 29634, USA Christopher S. McMahan††footnotemark:
Russell T. Shinohara Department of Biostatistics, Epidemiology, and
Informatics, Perelman School of Medicine, University of Pennsylvania,
Philadelphia, PA 19104, USA Kristin A. Linn††footnotemark:
for the Alzheimer’s Disease Neuroimaging Initiative111Data used in preparation
of this article were obtained from the Alzheimer’s Disease Neuroimaging
Initiative (ADNI) database (adni.loni.usc.edu). As such, the investigators
within the ADNI contributed to the design and implementation of ADNI and/or
provided data but did not participate in analysis or writing of this report. A
complete listing of ADNI investigators can be found at:
http://adni.loni.usc.edu/wp-
content/uploads/how_to_apply/ADNI_Acknowledgement_List.pdf

###### Abstract

Many analyses of neuroimaging data involve studying one or more regions of
interest (ROIs) in a brain image. In order to do so, each ROI must first be
identified. Since every brain is unique, the location, size, and shape of each
ROI varies across subjects. Thus, each ROI in a brain image must either be
manually identified or (semi-) automatically delineated, a task referred to as
segmentation. Automatic segmentation often involves mapping a previously
manually segmented image to a new brain image and propagating the labels to
obtain an estimate of where each ROI is located in the new image. A more
recent approach to this problem is to propagate labels from multiple manually
segmented atlases and combine the results using a process known as label
fusion. To date, most label fusion algorithms either employ voting procedures
or impose prior structure and subsequently find the maximum a posteriori
estimator (i.e., the posterior mode) through optimization. We propose using a
fully Bayesian spatial regression model for label fusion that facilitates
direct incorporation of covariate information while making accessible the
entire posterior distribution. We discuss the implementation of our model via
Markov chain Monte Carlo and illustrate the procedure through both simulation
and application to segmentation of the hippocampus, an anatomical structure
known to be associated with Alzheimer’s disease.

Key Words: Alzheimer’s disease, chromatic Gibbs sampling, conditionally
autoregressive model, image segmentation, latent variable modeling, regions of
interest

##  1 Introduction

The past twenty years have seen a dramatic increase in research involving
medical imaging, thanks in part to the advent of noninvasive techniques such
as magnetic resonance imaging (MRI). Neuroimaging in particular has benefited
from MRI technology due to the ability to delineate brain structures and
pathology with high spatial resolution. MRI technology itself is evolving and
improving, with new protocols and higher field strengths that lead to finer
resolution images with better contrast. With improvements to MRI technology
comes the challenge to identify and study smaller subregions of the brain,
known as regions of interest (ROI). Identification of ROIs is fundamental to
neuroimaging data analysis, as many studies involve analyzing population-level
structural variation and changes associated with disease, as well as
associations between ROI characteristics and clinical outcomes.

The manual delineation of medical images into ROIs by a trained expert is
often time-intensive and costly. Thus, even for moderate-sized research
studies, obtaining manual segmentations for each image can be infeasible
(Iglesias and Sabuncu, 2015). These so-called “gold standard” segmentations
provided by human experts, and the resulting ROI volumes, are known to exhibit
intra- and inter-rater variability (Yushkevich et al., 2006). In light of
these limitations, methods for automatic and semi-automatic segmentation of
ROIs are critical for increasing the feasibility of imaging studies.

A voxel is the single observational unit in a three-dimensional image. Each
voxel in an image contains a single grayscale intensity, so that an image can
be represented as an array of numeric values. An atlas is an image for which a
labeling exists that associates to each voxel a single label indicating the
structure or region to which it belongs. Early atlas-guided segmentation
techniques generally employed a single atlas to segment a single target or set
of target images (Pham et al., 2000). If one could obtain a perfectly
segmented atlas and a perfect registration (mapping) of the atlas to a target
image, a single-atlas approach would be sufficient. In practice, the
anatomical variability of even healthy brains cannot be captured by using a
single atlas. As a result, methods for using information from multiple atlases
have been introduced, including best atlas selection (Rohlfing et al., 2004)
and mapping images to a single atlas of label probabilities (Ashburner and
Friston, 2005). These approaches gave way to multi-atlas segmentation (MAS),
an approach that has been successfully employed in many application areas.
Doshi et al. (2016) propose MAS via ensembles of different registration
algorithms, as opposed to applying the same registration algorithm to every
atlas. Huo et al. (2016) propose the MaCRUISE algorithm that simultaneously
performs MAS and cortical surface reconstruction. Applications of MAS include
whole brain parcellation (Wang et al., 2012) and delineating individual
cortical/subcortical ROIs such as the cerebellum (Park et al., 2014), amygdala
(Hanson et al., 2012), and the ventricles (Raamana et al., 2014). A
particularly important application is segmentation and volumetry of the
hippocampus, atrophy of which is known to be associated with dementia and
Alzheimer’s disease (Yushkevich et al., 2006; Iglesias et al., 2013). Iglesias
and Sabuncu (2015) review MAS methods. Recently, deep learning also has been
applied to accomplish brain image segmentation. Examples include Akkus et al.
(2017), Milletari et al. (2017), and Wachinger et al. (2018).

MAS is characterized by combining information from a collection of atlases for
segmenting a novel, unlabeled image. The main components of MAS are 1) atlas
generation, 2) image registration and label propagation, and 3) label fusion.
Developments in the literature have largely focused on improving one or more
of these three components. Improvements in segmentation accuracy over single-
atlas methods are attributable to the wider range of anatomical variation that
is represented and leveraged by MAS. Our work here is focused on the label
fusion portion of MAS and assumes that multiple candidate segmentations are
available in the target image space. For example, these could have been
obtained by registering multiple manually segmented atlases to the target
image via nonlinear transformation and propagating each atlas’ labels to the
target image space.

A popular approach to label fusion is majority voting (Klein et al., 2005;
Heckemann et al., 2006), which takes the mode of the label distribution at
each voxel after registration. An early extension to simple majority voting is
weighted voting, where each atlas’ ‘vote’ is assigned a weight based on its
estimated reliability. These weights may be determined by mutual information
(Artaechevarria et al., 2008) or iterative re-weighting based on estimated
performance of each atlas (Langerak et al., 2010). Local weighting was
introduced to account for the fact that the quality of a given registration
can vary spatially over the image. Common approaches include weighting by the
voxelwise absolute image intensity differences (Isgum et al., 2009) or the
Jacobian of the registration transformation (Ramus and Malandain, 2010). Joint
label fusion (JLF; Wang et al., 2013) is currently considered the state-of-
the-art MAS method. JLF finds MAP estimates of atlas weights by exploiting the
fact that different atlases may be correlated in terms of their label errors.

Bayesian approaches to label fusion have also been proposed in the literature.
One of the first is that proposed by Sabuncu et al. (2010), which posits a
latent membership field to which each voxel belongs. A similar approach is
taken by the STAPLE algorithm (Warfield et al., 2004) in which the observed
rater segmentations are treated as corrupted versions of an underlying true
segmentation. Akhondi-Asl et al. (2014) extend the STAPLE algorithm by
associating estimated reliabilities with the membership field.

To date, most Bayesian approaches to label fusion are treated as optimization
problems to find the maximum a posteriori (MAP) estimator (i.e., the posterior
mode). They find the mode as a point estimate but stop short of quantifying
the uncertainty associated with the estimated classification. The MAP estimate
can be misleading, though, especially when the posterior distribution is
widely dispersed or strongly skewed (Fox and Nicholls, 2001). It is often the
case that nuisance parameters are estimated and held fixed instead of being
integrated out, thus not fully accounting for all sources of uncertainty in
the subsequent solution. Further, it is desirable to access the entire
posterior distribution to not only estimate a posterior probability map for
segmentation, but to have available any desired summary measures along with
appropriate measures of variability.

In this work, we propose a fully Bayesian model for binary label fusion of MR
images. We construct a latent variable regression model that incorporates
covariate information. We accommodate the variability of image registration
quality both within and between rater images by using spatially-varying
sensitivity and specificity processes, thereby facilitating smooth, local
weighting of each image as an inherent part of the model. We discuss prior
elicitation and implementation of the model via Markov chain Monte Carlo.

We motivate and present our proposed approach in Section 2, along with
practical considerations such as possible covariate information, prior
elicitation, and implementation. In Section 3 we consider a simulated example
with corrupted segmentations to compare its performance to that of simple and
weighted majority voting. We consider in Section 4 clinical MRI data obtained
from both healthy subjects and those that have been diagnosed with Alzheimer’s
disease. There we compare our proposed approach to voting procedures and JLF,
and illustrate how our approach can be used to estimate the distribution of
plausible volumes of a hippocampus. Section 5 concludes with a discussion of
the results and thoughts about future research directions.

##  2 Methods

###  2.1 Model

Suppose we are interested in segmenting a new brain image, hereinafter
referred to as the target image. We have available $R$ previously manually
segmented images to use as atlases after aligning the labeled brains to the
target. This alignment is done via image registration, the process of
transforming one image to the space of another via an affine or nonlinear map.
The registration problem is itself an ongoing area of research in the imaging
literature, but state of the art tools include the FLIRT and FNIRT functions
in FSL (Jenkinson and Smith, 2001; Jenkinson et al., 2002; Greve and Fischl,
2009; Jenkinson et al., 2012), DRAMMS (Ou et al., 2011), LDDMM (e.g., Zhong et
al., 2010), and SyN registration (Tustison and Avants, 2013). In what follows,
we assume that a single algorithm is used to register each atlas to the
target, but our proposed framework could easily account for different
registration algorithms. In either case, the sources of uncertainty leading to
corruption of the atlas labels are the radiologists’ imperfect manual
segmentations, the set of brains that have been manually segmented, none of
which is identical to the target brain, and the registration algorithm itself,
which cannot produce a perfect alignment between images.

After registering each atlas to the target, we have $R$ images indexed by
$v=1,\ldots,V$, where $V$ is the number of voxels. Corresponding to atlas
$r\in\\{1,\ldots,R\\}$ is a set of labels
$\\{Y_{1r},\ldots,Y_{Vr}\\}\in\\{0,1\\}^{V}$. The atlases will disagree with
each other after registration and will each be a corrupted version of the
underlying truth. In addition to the observed labels $Y_{vr}$, let
$T_{v}\in\\{0,1\\},~{}v=1,\ldots,V$, denote the “true” (unobservable) status
of the voxel $v$, so that $T_{v}=1$ indicates that voxel $v$ is a part of the
structure of interest, and $T_{v}=0$ otherwise.

We associate to each atlas $r$ a sensitivity $\xi$, an ability to correctly
detect when a voxel is truly part of the structure of interest (i.e., to avoid
false negatives). Similarly, we have also for each atlas a specificity $\psi$,
an ability to correctly determine when a voxel is truly outside of the ROI
(i.e., to avoid a false positive). Since the quality of any registration can
vary throughout an image, the best atlas to use depends on location. In other
words, poor agreement between the atlas and a target image in one part of the
brain does not necessarily mean that the atlas is unreliable in other areas
(Artaechevarria et al., 2009). Thus, we allow the sensitivity and specificity
of each atlas to be a function not only of the specific atlas, but also the
voxel. We suppose that $P(Y_{vr}=1\mid T_{v}=1)=\xi(v,r)$ and $P(Y_{vr}=0\mid
T_{v}=0)=\psi(v,r)$, for atlases $r=1,\ldots,R$ and voxels $v=1,\ldots,V$.
Given the true voxel statuses and the sensitivity and specificity of each
atlas, we assume that observed labels are conditionally independent,

$Y_{vr}\mid T_{v},~{}\xi(v,r),~{}\psi(v,r)\stackrel{{\scriptstyle
ind.}}{{\sim}}\text{Bern}(p^{\ast}(v,r)),~{}~{}v=1,\ldots,V;~{}r=1,\ldots R,$
(1)

where $p^{\ast}(v,r)=\xi(v,r)^{T_{v}}[1-\psi(v,r)]^{1-T_{v}}$. In particular,
our model assumes that any dependence between $Y_{vr}$ and
$Y_{v^{\prime}r^{\prime}}$, for $v\neq v^{\prime}$ or $r\neq r^{\prime}$, is
completely explained by the reliability of the atlases. For example, if
$Y_{vr_{0}}$ and $Y_{v^{\prime}r_{0}}$ come from the same atlas $r_{0}$, we
assume the dependence between them is only due to the dependence between
$\xi(v,r_{0})$ and $\xi(v^{\prime},r_{0})$ and between $\psi(v,r_{0})$ and
$\psi(v^{\prime},r_{0})$.

We propose a framework for incorporating spatial smoothness into the
sensitivity and specificity processes, as well as additional covariate
information that can be informative with respect to the true voxel status
$T_{v}$. The spatial process and covariate models themselves depend on
unobservable and unkonwn parameters. Hence we take a hierarchical Bayesian
approach and assign them prior distributions as well. The prior distributions
and concomitant hyper-parameters are described in this Subsection.

It is important to impose smoothness on the reliabilities to mitigate the
effect of image noise (Artaechevarria et al., 2009). We model the
sensitivities and specificities as

$\displaystyle\begin{aligned}
\xi(v,r)&=\Phi(\boldsymbol{x}_{vr}^{\top}\boldsymbol{\beta}_{r}+\phi_{vr})\\\
\psi(v,r)&=\Phi(\boldsymbol{z}_{vr}^{\top}\boldsymbol{\gamma}_{r}+\eta_{vr}),\end{aligned}$
(2)

where $\Phi(\cdot)$ is the standard Gaussian CDF,
$\boldsymbol{x}_{vr}\in\mathbb{R}^{\ell}$ and
$\boldsymbol{z}_{vr}\in\mathbb{R}^{k}$ are known vectors of covariates with
coefficients $\boldsymbol{\beta}_{r}$ and $\boldsymbol{\gamma}_{r}$,
respectively, and $\phi_{vr}$ and $\eta_{vr}$ are elements of mean-zero
Gaussian Markov random fields (GMRFs; Rue and Held, 2005). Our experience is
that it is often sufficient to simply set
$\boldsymbol{x}=\boldsymbol{z}=\boldsymbol{0}$ to let the spatial processes
detect the trends a posteriori. Observe that we assume that $\xi(v,r)$ and
$\psi(v,r)$ are independent of each other, since they are defined by two
mutually exclusive outcomes, $T_{v}=1$ and $T_{v}=0$.

We assume
$\boldsymbol{\phi}_{r}=(\phi_{1r},...,\phi_{Vr})^{\top}\in\mathbb{R}^{V}$ and
$\boldsymbol{\eta}_{r}=(\eta_{1r},...,\eta_{Vr})^{\top}\in\mathbb{R}^{V}$,
$r=1,...,R$, each independently follow conditionally autoregressive models
(CAR; Besag, 1974), given by

$\displaystyle\begin{aligned} \boldsymbol{\phi}_{r}&\stackrel{{\scriptstyle
ind.}}{{\sim}}N_{V}\left\\{\boldsymbol{0},\tau_{\boldsymbol{\phi}_{r}}^{-1}(\boldsymbol{D}-\rho_{\boldsymbol{\phi}_{r}}\boldsymbol{W})^{-1}\right\\}\\\
\boldsymbol{\eta}_{r}&\stackrel{{\scriptstyle
ind.}}{{\sim}}N_{V}\left\\{\boldsymbol{0},\tau_{\boldsymbol{\eta}_{r}}^{-1}(\boldsymbol{D}-\rho_{\boldsymbol{\eta}_{r}}\boldsymbol{W})^{-1}\right\\}.\end{aligned}$
(3)

Here, $\boldsymbol{W}=\\{w_{ij}\\}_{i,j=1}^{V}\in\mathbb{R}^{V\times V}$ is a
neighborhood matrix such that $w_{v,v^{\prime}}=I(v\sim v^{\prime})$, where
$v\sim v^{\prime}$ if and only if voxels $v$ and $v^{\prime}$ share an edge or
a corner, $I(\cdot)$ is the indicator function, and
$\boldsymbol{D}=\text{diag}\left(\sum_{j=1}^{V}w_{ij},~{}i=1,\ldots,V\right)\in\mathbb{R}^{V\times
V}$. This yields a nonstationary process due to the edge effects; i.e., the
voxels on the edges have different numbers of neighbors than the interior
voxels, leading to different marginal variances. Besag and Kooperberg (1995)
observe that when images are of large dimension and the regions of interest
are in the interior, edge effects are negligible with respect to inferences,
so one can safely ignore them.

When $T_{v}=0$ for all $v$, $\xi(v,r)$ does not appear in the likelihood
determined by (1) and hence $\boldsymbol{\phi}_{r}$ is not updated in the
posterior. Since $T_{1}=\cdots=T_{V}=0$ with positive probability, the prior
distribution on $\boldsymbol{\phi}_{r}$ must be proper to avoid an improper
posterior, and likewise for $\boldsymbol{\eta}_{r}$. Thus, we include the
“propriety parameters” $\rho_{\boldsymbol{\phi}_{r}}$ and
$\rho_{\boldsymbol{\eta}_{r}}$ in (3) to force the precision matrices to be
nonsingular, as suggested by Banerjee et al. (2015). A sufficient condition is
$|\rho_{\boldsymbol{\phi}_{r}}|<1$, in which case
$(\boldsymbol{D}-\rho_{\boldsymbol{\phi}_{r}}\boldsymbol{W})$ is diagonally
dominant and hence positive definite.

It is likely that we will have available auxiliary information concerning the
true location of the structure of interest. Toward this end, we suppose that,

$P(T_{v}=1\mid\boldsymbol{\delta})=g^{-1}(\boldsymbol{c}_{v}^{\top}\boldsymbol{\delta}),~{}~{}v=1,\ldots,V,$
(4)

where $\boldsymbol{c}_{v}\in\mathbb{R}^{J}$ is a vector of covariates
pertaining to voxel $v$, $\boldsymbol{\delta}$ is the corresponding vector of
regression coefficients, and $g:(0,1)\rightarrow\mathbb{R}$ is a one-to-one
and differentiable link function. We assume any dependence between inclusion
indicators is completely explained by the covariate information; i.e.,
$T_{v}\perp T_{v^{\prime}}\mid\boldsymbol{\delta},~{}v\neq v^{\prime}$. In
Sections 3 and 4, we use signed distance label maps and either normalized
image intensity or tissue class segmentations as covariate information for
$\boldsymbol{c}_{v}$. In Subsection 2.2, we discuss an approach for eliciting
a prior distribution on $\boldsymbol{\delta}$.

The model is completed with prior distributions on the regression coefficients
and precision parameters in (2), (3), and (4). We take conventional Gaussian
and Gamma priors for these parameters, since they are sufficiently flexible
with appropriate specification of the corresponding hyperparameters. Thus, we
assume

$\displaystyle\begin{aligned} \boldsymbol{\beta}_{r}\stackrel{{\scriptstyle
iid}}{{\sim}}N_{\ell}\left(\boldsymbol{0},\boldsymbol{\Sigma}_{\boldsymbol{\beta}}\right);~{}~{}\boldsymbol{\gamma}_{r}&\stackrel{{\scriptstyle
iid}}{{\sim}}N_{k}\left(\boldsymbol{0},\boldsymbol{\Sigma}_{\boldsymbol{\gamma}}\right),~{}r=1,\ldots,R\\\
\boldsymbol{\delta}&{\sim}N_{J}\left(\boldsymbol{0},\boldsymbol{\Sigma}_{\boldsymbol{\delta}}\right),\\\
\tau_{\boldsymbol{\phi}_{r}}\stackrel{{\scriptstyle
iid}}{{\sim}}Ga\left(a_{\boldsymbol{\phi}},b_{\boldsymbol{\phi}}\right);~{}~{}\tau_{\boldsymbol{\eta}_{r}}&\stackrel{{\scriptstyle
iid}}{{\sim}}Ga\left(a_{\boldsymbol{\eta}},b_{\boldsymbol{\eta}}\right),~{}r=1,\ldots,R.\end{aligned}$
(5)

Let $\boldsymbol{\omega}$ denote the collection of all model parameters except
the inclusion indicators
$\boldsymbol{T}=(T_{1},\ldots,T_{V})^{\top}\in\mathbb{R}^{V}$. Under the model
determined by equations (1) - (5), the joint posterior density of the
parameters, conditional on the observed data
$\boldsymbol{Y}=(Y_{11},\ldots,Y_{VR})^{\top}\in\mathbb{R}^{VR}$, is

$\displaystyle\pi(\boldsymbol{T},\boldsymbol{\omega}\mid\boldsymbol{Y})$
$\displaystyle\propto$
$\displaystyle\prod_{v=1}^{V}\prod_{r=1}^{R}\Phi(\phi_{vr})^{Y_{vr}T_{v}}\\{1-\Phi(\phi_{vr})\\}^{(1-Y_{vr})T_{v}}\Phi(\eta_{vr})^{(1-Y_{vr})(1-T_{v})}\\{1-\Phi(\eta_{vr})\\}^{Y_{vr}(1-T_{v})}$
$\displaystyle\times\prod_{v=1}^{V}\\{g^{-1}(\boldsymbol{c}_{v}^{\top}\boldsymbol{\delta})\\}^{T_{v}}\\{1-g^{-1}(\boldsymbol{c}_{v}^{\top}\boldsymbol{\delta})\\}^{1-T_{v}}\exp\left(-\boldsymbol{\delta}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\delta}}^{-1}\boldsymbol{\delta}/2\right)$
$\displaystyle\times\prod_{r=1}^{R}\tau_{\boldsymbol{\phi}_{r}}^{V/2}\exp\left\\{-\tau_{\boldsymbol{\phi}_{r}}(\boldsymbol{\phi}_{r}-\boldsymbol{X}_{r}\boldsymbol{\beta}_{r})^{\top}(\boldsymbol{D}-\rho_{\boldsymbol{\phi}_{r}}\boldsymbol{W})(\boldsymbol{\phi}_{r}-\boldsymbol{X}_{r}\boldsymbol{\beta}_{r})/2\right\\}$
$\displaystyle\times\prod_{r=1}^{R}\tau_{\boldsymbol{\eta}_{r}}^{V/2}\exp\left\\{-\tau_{\boldsymbol{\eta}_{r}}(\boldsymbol{\eta}_{r}-\boldsymbol{Z}_{r}\boldsymbol{\gamma}_{r})^{\top}(\boldsymbol{D}-\rho_{\boldsymbol{\eta}_{r}}\boldsymbol{W})(\boldsymbol{\eta}_{r}-\boldsymbol{Z}_{r}\boldsymbol{\gamma}_{r})/2\right\\}$
$\displaystyle\times\prod_{r=1}^{R}\textrm{exp}\left(-\boldsymbol{\beta}_{r}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\beta}}^{-1}\boldsymbol{\beta}_{r}/2\right)\exp\left(-\boldsymbol{\gamma}_{r}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\gamma}}^{-1}\boldsymbol{\gamma}_{r}/2\right)$
$\displaystyle\times\prod_{r=1}^{R}\tau_{\boldsymbol{\phi}_{r}}^{a_{\boldsymbol{\phi}_{r}}-1}\tau_{\boldsymbol{\eta}_{r}}^{a_{\boldsymbol{\eta}_{r}}-1}\exp(-\tau_{\boldsymbol{\phi}_{r}}b_{\boldsymbol{\phi}_{r}})\exp(-\tau_{\boldsymbol{\eta}_{r}}b_{\boldsymbol{\eta}_{r}}).$

To aid with exposition, Figure 1 provides a graphical depiction of our
proposed model.

Figure 1: Graphical depiction of the proposed Bayesian label fusion model for
segmenting brain images consisting of $V$ voxels using $R$ available atlases.
The observed data are $Y_{vr}$, the latent true inclusion indicators are
denoted $T_{v}$, and for each rater $r$ and voxel $v$, there is a sensitivity
$\xi(v,r)$ and specificity $\psi(v,r)$. Note that CMP indicates the
conditional mean prior specification (Bedrick et al., 1996), and GMRF
indicates a Gaussian Markov random field.

###  2.2 Prior Specification

To specify the hyperparameters determining the prior distributions for
$\boldsymbol{\beta}_{r},\boldsymbol{\gamma}_{r}$, $\boldsymbol{\phi}_{r}$, and
$\boldsymbol{\eta}_{r}$ appearing in (2), we recognize that
$\Phi(4)-\Phi(-4)\approx 1$, so that the effect of the mean function on the
likely values of the reliability parameters is almost certainly between $-4$
and $4$. For instance, considering the sensitivity, a priori we suppose that
for any $v$ and any $r$,
$P(|\boldsymbol{x}_{vr}^{\top}\boldsymbol{\beta}_{r}+\phi_{vr}|\leq 4)\approx
1$ for all $\boldsymbol{x}_{vr}$. With a pure spatial process
($\boldsymbol{x}_{vr}=\boldsymbol{0}$, for all $v,r$), this requirement
becomes $P(|\phi_{vr}|\leq 4)\approx 1$. To use this to induce a prior on the
spatial effect $\phi_{vr}$, we use the fact that
$Var(\phi_{vr})\approx(0.7^{2}w_{v}.\tau_{\boldsymbol{\phi}_{r}})^{-1}$, where
$w_{v}.=\sum_{k=1}^{V}w_{vk}$ (Bernardinelli et al., 1995; Eberly and Carlin,
2000). With eight neighbors at a voxel $v$, the desired constraint can be
solved to yield $\tau_{\boldsymbol{\phi}_{r}}\approx 0.5$. Hence we take
$a_{\phi}=1$ and $b_{\phi}=2$ in (5). The argument is the same for specifying
the hyperparameters on $\boldsymbol{\gamma}_{r}$ and
$\tau_{\boldsymbol{\eta}_{r}}$.

Similar to specifying the priors on the reliability parameters, we take
$\boldsymbol{\Sigma}_{\boldsymbol{\delta}}=d\boldsymbol{I}$ in (5), for some
appropriate $d>0$ and $\boldsymbol{I}\in\mathbb{R}^{J\times J}$ the identity
matrix, to allow a reasonably flat prior on the $p_{v}:=P(T_{v}=1)$ scale.
Alternatively, with prior knowledge concerning the underlying structure of
interest, it is possible to induce a prior on $\boldsymbol{\delta}$ through a
so-called conditional mean prior (CMP; Bedrick et al., 1996). Under this
approach, we partition the $J$-dimensional covariate space into $J$ regions
and choose covariate vectors
$\widetilde{\boldsymbol{c}}_{1},\widetilde{\boldsymbol{c}}_{2},\ldots,\widetilde{\boldsymbol{c}}_{J}$
so that
$\widetilde{\boldsymbol{C}}=(\widetilde{\boldsymbol{c}}_{1}~{}\widetilde{\boldsymbol{c}}_{2}~{}\ldots\widetilde{\boldsymbol{c}}_{J})^{\top}\in\mathbb{R}^{J\times
J}$ is nonsingular. A prior is put directly on the mean response at these
covariate values,
$E[(\widetilde{T}_{1},\ldots,\widetilde{T}_{J})^{\top}\mid\boldsymbol{\delta}]=\boldsymbol{G}^{\prime}(\widetilde{\boldsymbol{C}}\boldsymbol{\delta})~{}\in\mathbb{R}^{J}$,
where $\boldsymbol{G}^{\prime}(\cdot)$ applies $g^{-1}$ componentwise. While
the regression coefficients themselves may be difficult to interpret, we can
meaningfully assign a distribution to
$\widetilde{p}_{j}=g^{-1}(\widetilde{\boldsymbol{c}}_{j}^{\top}\boldsymbol{\delta})=P(\widetilde{T}_{j}=1\mid\boldsymbol{\delta})$
using a Beta distribution with shape parameters chosen to reflect the prior
knowledge at the covariate values. Bedrick et al. (1996) show that if
$\widetilde{p}_{j}\stackrel{{\scriptstyle
indep.}}{{\sim}}\textrm{Beta}(a_{\widetilde{p}_{j}},b_{\widetilde{p}_{j}}),~{}j=1,\ldots,J$,
then the induced prior is
$\pi(\boldsymbol{\delta})\propto\prod_{j=1}^{J}\\{g^{-1}(\widetilde{\mathbf{c}}_{j}^{\top}\boldsymbol{\delta})\\}^{a_{\widetilde{p}_{j}}-1}\\{1-g^{-1}(\widetilde{\mathbf{c}}_{j}^{\top}\boldsymbol{\delta})\\}^{b_{\widetilde{p}_{j}}-1}\dot{g}^{-1}(\widetilde{\mathbf{c}}_{j}^{\top}\boldsymbol{\delta})$.

###  2.3 Implementation via MCMC

For posterior inference, we use Markov chain Monte Carlo sampling (MCMC;
Gelfand and Smith, 1990), specifically a Metropolis-Hastings-within-Gibbs
scheme (Metropolis et al., 1953; Hastings, 1970; Geman and Geman, 1984;
Müller, 1991). The necessary full conditional distributions are given in
Section 1 of the Supplementary Material.

Sampling $\boldsymbol{\delta}$ can be done using data augmentation (Albert and
Chib, 1993). An alternative is the algorithm of Gamerman (1997). Under our
proposed model (1) - (5), the full conditional density of
$\boldsymbol{\delta}$ is
$\pi(\boldsymbol{\delta}|\boldsymbol{T},\boldsymbol{C})\propto\exp\\{-\boldsymbol{\delta}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\delta}}^{-1}\boldsymbol{\delta}/2+$
$\sum_{v=1}^{V}T_{v}\theta_{v}-b(\theta_{v})\\}$, where
$\boldsymbol{C}\in\mathbb{R}^{V\times J}$ is the design matrix constructed
from $\boldsymbol{c}_{1},\ldots,\boldsymbol{c}_{V}$,
$\theta_{v}=\log[g^{-1}(\boldsymbol{c}_{v}^{\top}\boldsymbol{\delta})/\\{1-g^{-1}(\boldsymbol{c}_{v}^{\top}\boldsymbol{\delta})\\}],$
and $b(\cdot)=\log\\{1+\exp(\cdot)\\}$. Define
$\boldsymbol{T}(\boldsymbol{\delta})=(T_{1}(\boldsymbol{\delta}),\ldots,$
$T_{V}(\boldsymbol{\delta}))^{\top}\in\mathbb{R}^{V}$ and
$\boldsymbol{Q}(\boldsymbol{\delta})=\text{diag}(Q_{1}(\boldsymbol{\delta}),\ldots,Q_{V}(\boldsymbol{\delta}))^{\top}$
with
$T_{v}(\boldsymbol{\delta})=\boldsymbol{c}_{v}^{\top}\boldsymbol{\delta}+\\{T_{v}-g^{-1}(\mathbf{c}_{v}^{\top}\boldsymbol{\delta})\\}\dot{g}(g^{-1}(\boldsymbol{c}_{v}^{\top}\boldsymbol{\delta}))$
and
$Q_{v}^{-1}(\boldsymbol{\delta})=\ddot{b}(\theta_{v})\\{\dot{g}(g^{-1}(\boldsymbol{c}_{v}^{\top}\boldsymbol{\delta}))\\}^{2}$
for $v=1,\ldots,V$. Then, given the current iterate
$\boldsymbol{\delta}^{(t)}$, the proposal distribution is
$\boldsymbol{\delta}^{(t+1)}\mid\boldsymbol{\delta}^{(t)}\sim
N_{{\color[rgb]{0,0,1}J}}(\boldsymbol{m}^{(t)},\boldsymbol{V}^{(t)}),$ where
$\boldsymbol{m}^{(t)}=\\{\boldsymbol{\Sigma}_{\boldsymbol{\delta}}^{-1}+\boldsymbol{C}^{\top}\boldsymbol{Q}(\boldsymbol{\delta}^{(t))})\boldsymbol{C}\\}^{-1}$
$\boldsymbol{C}^{\top}\boldsymbol{Q}(\boldsymbol{\delta}^{(t)})\boldsymbol{T}(\boldsymbol{\delta}^{(t)})$
and
$\boldsymbol{V}^{(t)}=\\{\boldsymbol{\Sigma}_{\boldsymbol{\delta}}^{-1}+\boldsymbol{C}^{\top}\boldsymbol{Q}(\boldsymbol{\delta}^{(t)})\boldsymbol{C}\\}^{-1}$.
If the prior on $\boldsymbol{\delta}$ is elicited through the CMP approach,
this proposal mechanism can still be used after suitable modification. See
Section 2 of the Supplementary Material for a derivation under the logistic
link.

Depending on $\boldsymbol{T}$, every voxel contributes information to either
$\boldsymbol{\phi}_{r}$ or $\boldsymbol{\eta}_{r}$, but not both, due to the
definitions of $\xi(\cdot,r)$ and $\psi(\cdot,r)$. Following Albert and Chib
(1993), we introduce latent variables $\zeta_{vr}^{1}\mid
T_{v},\boldsymbol{\beta}_{r},\phi_{vr}\stackrel{{\scriptstyle
ind}}{{\sim}}N(T_{v}(\boldsymbol{x}_{vr}^{\top}\boldsymbol{\beta}_{r}+\phi_{vr}),1),~{}\zeta_{vr}^{0}\mid
T_{v},\boldsymbol{\gamma}_{r},\eta_{vr}\stackrel{{\scriptstyle
ind}}{{\sim}}N((1-T_{v})(\boldsymbol{z}_{vr}^{\top}\boldsymbol{\gamma}_{r}+\eta_{vr}),1),~{}v=1,\ldots,V$.
Letting $\mathcal{V}_{\ell}=\\{v:T_{v}=\ell\\},~{}\ell=0,1,$ define
$Y_{vr}=I(\zeta_{vr}^{1}\geq 0)$ for $v\in\mathcal{V}_{1}$, and similarly
$Y_{vr}=I(\zeta_{vr}^{0}<0)$ for $v\in\mathcal{V}_{0}$. Letting
$\boldsymbol{\mathcal{T}}=\text{diag}(T_{1},\ldots,T_{V})~{}\in\mathbb{R}^{V\times
V}$,
$\boldsymbol{\zeta}_{r}^{1}=(\zeta_{1r}^{1},\ldots,\zeta_{Vr}^{1})^{\top}\in\mathbb{R}^{V}$,
and
$\boldsymbol{\zeta}_{r}^{0}=(\zeta_{1r}^{0},\ldots,\zeta_{Vr}^{0})^{\top}\in\mathbb{R}^{V}$,
this can be expressed as
$\boldsymbol{\zeta}_{r}^{1}\mid\boldsymbol{\beta}_{r},\phi_{vr},\boldsymbol{T}\sim
N_{V}\\{\boldsymbol{\mathcal{T}}(\boldsymbol{\phi}_{r}+\boldsymbol{X}_{r}\boldsymbol{\beta}_{r}),\boldsymbol{I}\\}$
and
$\boldsymbol{\zeta}_{r}^{0}\mid\boldsymbol{\beta}_{r},\phi_{vr},\boldsymbol{T}\sim
N_{V}\\{(\boldsymbol{I}-\boldsymbol{\mathcal{T}})(\boldsymbol{\eta}_{r}+\boldsymbol{Z}_{r}\boldsymbol{\gamma}_{r}),\boldsymbol{I}\\}$.
Using the augmented data likelihood, the full conditional distributions of the
spatial effects are given by
$\boldsymbol{\phi}_{r}|\boldsymbol{\zeta}_{r}^{1},\tau_{\boldsymbol{\phi}_{r}},\boldsymbol{\beta}_{r}\sim
N_{{\color[rgb]{0,0,1}V}}(\boldsymbol{\mu}_{\boldsymbol{\phi}_{r}}^{*},\boldsymbol{\Sigma}_{\boldsymbol{\phi}_{r}}^{*})$
and
$\boldsymbol{\eta}_{r}|\boldsymbol{\zeta}_{r}^{0},\tau_{\boldsymbol{\eta}_{r}},\boldsymbol{\gamma}_{r}\sim
N_{{\color[rgb]{0,0,1}V}}(\boldsymbol{\mu}_{\boldsymbol{\eta}_{r}}^{*},\boldsymbol{\Sigma}_{\boldsymbol{\eta}_{r}}^{*})$,
for $r=1,...,R$, where
$\boldsymbol{\Sigma}_{\boldsymbol{\phi}_{r}}^{*}=\left(\boldsymbol{\mathcal{T}}+\tau_{\boldsymbol{\phi}_{r}}(\boldsymbol{D}-\rho_{\phi}\boldsymbol{W})\right)^{-1},$
$\boldsymbol{\Sigma}_{\boldsymbol{\eta}_{r}}^{*}=\left(\boldsymbol{I}-\boldsymbol{\mathcal{T}}+\tau_{\boldsymbol{\eta}_{r}}(\boldsymbol{D}-\rho_{\eta}\boldsymbol{W})\right)^{-1}$,
$\boldsymbol{\mu}_{\boldsymbol{\phi}_{r}}^{*}=\boldsymbol{\Sigma}_{\boldsymbol{\phi}_{r}}^{*}\boldsymbol{\mathcal{T}}^{\top}(\boldsymbol{\zeta}_{r}^{1}-\boldsymbol{\mathcal{T}}\boldsymbol{X}_{r}\boldsymbol{\beta}_{r})$,
and
$\boldsymbol{\mu}_{\boldsymbol{\eta}_{r}}^{*}=\boldsymbol{\Sigma}_{\boldsymbol{\eta}_{r}}^{*}(\boldsymbol{I}-\boldsymbol{\mathcal{T}})^{\top}(\boldsymbol{\zeta}_{r}^{0}-(\boldsymbol{I}-\boldsymbol{\mathcal{T}})\boldsymbol{Z}_{r}\boldsymbol{\gamma}_{r})$.
For the conditional distribution of $\boldsymbol{\zeta}$, it can be shown
that, for all $v,r$, $\zeta_{vr}^{1}\mid
Y_{vr},\phi_{vr},\beta_{r},T_{v}\stackrel{{\scriptstyle
ind}}{{\sim}}I(Y_{vr}=0)TN\\{T_{v}(\phi_{vr}+\boldsymbol{x}_{vr}^{\top}\boldsymbol{\beta}_{r}),1,(-\infty,0)\\}+I(Y_{vr}=1)TN\\{T_{v}(\phi_{vr}+\boldsymbol{x}_{vr}^{\top}\boldsymbol{\beta}_{r}),1,(0,\infty)\\}$,
and $\zeta_{vr}^{0}\mid
Y_{vr},\eta_{vr},\beta_{r},T_{v}\stackrel{{\scriptstyle
ind}}{{\sim}}I(Y_{vr}=1)TN\\{(1-T_{v})(\eta_{vr}+\boldsymbol{z}_{vr}^{\top}\boldsymbol{\gamma}_{r}),1,(-\infty,0)\\}+I(Y_{vr}=0)TN\\{(1-T_{v})(\eta_{vr}+\boldsymbol{z}_{vr}^{\top}\boldsymbol{\gamma}_{r}),1,(0,\infty)\\}$,
where $TN(\mu,\sigma^{2},(a,b))$ denotes a Gaussian distribution with mean
$\mu$ and variance $\sigma^{2}$ truncated to the interval $(a,b)$.

A convenient feature of CAR models is the availability of easily interpretable
conditional distributions that follow from the Markov property. For instance,
the CAR prior on $\boldsymbol{\phi}_{r}$ implies that
$\phi_{vr}\mid\boldsymbol{\phi}_{(-v),r},\tau_{\phi_{r}}\sim
N(\rho\overline{\phi}_{v},(w_{v}.\tau_{\phi_{r}})^{-1})$, where
$\boldsymbol{\phi}_{(-v),r}=(\phi_{1r},\ldots,\phi_{v-1,r},\phi_{v+1,r},\ldots,$
$\phi_{Vr})^{\top}\in\mathbb{R}^{V-1}$, $w_{v}.=\sum_{k=1}^{V}w_{vk}$, and
$\overline{\phi}_{v}=w_{v\cdot}^{-1}\sum_{k=1}^{V}w_{vk}\phi_{kr}$. Combined
with the model for the augmented data $\zeta_{vr}^{1}$, this yields
$\phi_{vr}\mid\zeta_{vr}^{1},\boldsymbol{\phi}_{(-v),r},\tau_{\phi},\boldsymbol{\beta}_{r},T_{v}\sim
N(\mu_{v,r}^{\ast},\sigma_{v}^{2,\ast})$, where
$\sigma_{v}^{2,\ast}=(T_{v}+\tau_{\phi}w_{v}.)^{-1}$ and
$\mu_{v,r}^{\ast}=\sigma^{2.\ast}(T_{v}(\zeta_{vr}^{1}-\boldsymbol{x}_{vr}^{\top}\boldsymbol{\beta}_{r})+\rho_{\phi}w_{v}.\tau_{\phi}\overline{\phi}_{v})$.
The Markov property still holds since, conditional on its neighbors,
$\phi_{vr}$ is independent of the rest of the field. Thus, we can update
independent elements of $\boldsymbol{\phi}_{r}$ simultaneously (e.g., in
parallel or through ‘vectorized’ calculations in `R` (R Development Core Team,
2008)) while conditioning on all of their common neighbors. Such simultaneous
updating is accomplished through chromatic sampling (Gonzalez et al., 2011;
Brown et al., 2019), so named due to its association with coloring the graph
associated with the GMRF. With simultaneous updates, chromatic Gibbs sampling
dramatically reduces the cost of sequential single-site updating while
avoiding the computational burden of drawing from high-dimensional Gaussian
distributions. In fact, Brown et al. (2019) show that its computational cost
increases at a much lower rate than block sampling via multivariate Gaussian
draws. We use chromatic Gibbs sampling on both $\boldsymbol{\phi}_{r}$ and
$\boldsymbol{\eta}_{r}$ to accelerate the MCMC routines.

##  3 Numerical Studies

Here we illustrate our proposed approach on a simulated example. It may occur
in practice that multiple reference atlases are subject to the same
deleterious effect; e.g., they are all registered using a poor algorithm. This
could result in multiple atlases closely agreeing on an incorrect
segmentation. If we happen to have one atlas that registered well and provides
a good segmentation, it could still be outweighed by multiple poor atlases,
resulting in a poor segmentation when using a voting-based procedure.

To consider this possibility, we simulate a structure on a $40\times 40$ grid,
displayed in Figure 2. We simulate also image intensities as might be obtained
from an additional imaging modality. Also displayed in Figure 2 is one
reliable candidate segmentation and three poor segmentations, the latter of
which closely agree. Voxelwise image intensity similarity is sometimes used to
quantify rater reliability throughout an image (Isgum et al., 2009). As such,
we simulate image intensity similarities for each atlas, displayed in
Supplementary Figure 1. To make the assessment more challenging (and more
realistic), we slightly offset the intensity differences from the areas of
poor segmentations.

Figure 2: True structure of interest (far left), along with the associated
normalized image intensities (middle left), one quality segmentation (middle
right) and three poor atlas segmentations (far right) for the simulated label
fusion example. The three poor segmentations are superimposed to indicate
their agreement.

We find the signed distance label transforms for each atlas (Iglesias et al.,
2012) to use as a covariate. A signed distance label (SDL) transform for a
binary image assigns to each pixel a number corresponding to its distance from
the nearest edge in an image, where an edge is indicated by a zero adjacent to
a one. The sign of the distance corresponds to whether or not the pixel is
inside or outside an identified structure in the binary image; e.g., a pixel
has negative distance if it is inside the structure, positive distance if it
is outside, and zero if it is on the boundary. We take a weighted sum over the
atlas-specific SDL maps, where the weights are inversely proportional to the
square of the intensity difference between the rater image and the target;
i.e.,
$c_{v1}=\sum_{r=1}^{4}|i_{r}(v)-i_{t}(v)|^{-2}d_{r}(v)/(\sum_{r=1}^{4}|i_{r}(v)-i_{t}(v)|^{-2})$,
where $i_{r}(v)$ and $i_{t}(v)$ are the intensity values for rater image $r$
and the target image at pixel $v$, respectively, and $d_{r}(v)$ is the SDL
assigned to pixel $v$ from rater image $r$. We use the pixel-level intensities
as additional covariates, $c_{v2}=i_{t}(v)$. It is important to include a
distance-intensity interaction, $c_{v3}:=c_{v1}c_{v2}$, since high intensity
voxels far away from the object of interest should not be included, but the
information can be helpful when we believe a voxel is close to the object.

We use CMP to induce a prior on the $\boldsymbol{\delta}$ coefficients in (4).
(See Section 2.) The dimension of the predictor space is $J=4$, so we augment
the model with pseudo-observations under covariate values corresponding to (i)
small distance, high intensity; (ii) large distance, low intensity; (iii) mid-
distance and average intensity; and (iv) large distance but high intensity. We
impose prior knowledge that (i) is very likely to be included in the object,
(ii) is very unlikely to be included, (iii) is a borderline case, and (iv) is
not very likely to be included. (See Supplementary Material for more details.)
We take the GMRF precision parameters to be $a_{\phi}=a_{\eta}=1$ and
$b_{\phi}=b_{\eta}=2$.

For posterior inference, a single Monte Carlo Markov chain is run for 100,000
iterations, thinning to every $25^{\text{th}}$ iterate to save memory and
reduce autocorrelation. The first half of the chain is discarded as a burn-in
period. We use trace plots of the $\boldsymbol{\delta}$ coefficients along
with Geweke statistics (Geweke, 1992) and lag 1 autocorrelations as
convergence diagnostics. The diagnostics, displayed in Supplementary Figure 2,
provide evidence of convergence.

Figure 3 displays the posterior mean estimate of the target along with pixel-
wise standard deviations. To reduce the Monte Carlo variance, we use “Rao-
Blackwellized” estimators (Gelfand and Smith, 1990),
$\hat{P}(T_{v}=1\mid\boldsymbol{Y})=N^{-1}\sum_{k=1}^{N}\hat{P}(T_{v}=1\mid\boldsymbol{\omega}^{(k)},\boldsymbol{Y})$,
where $\boldsymbol{\omega}^{(k)}$ denotes the $k^{\text{th}}$ MCMC iterate of
all the parameters except $T_{1},\ldots,T_{V}$. Close agreement between the
target and the estimate is evident. As expected, the uncertainty about the
structure is largest near the edges of the observed segmentations, as well as
in the area of higher intensity. The effect of high intensity is strongly
attenuated when it is thought to be well away from the target object.
Posterior inference also produces estimates of the reliabilities of the rater
atlases. Supplementary Figure 3 displays the estimated sensitivity and
specificity fields for each of the four atlas segmentations. We see that all
atlases have high specificity in regions well away from the structure. Near
the object, however, we observe high sensitivity of the good atlas with lower
values for the poor ones. Conversely, the tendencies for the bad atlases to
produce false positives are seen in the areas of low specificity.

Figure 3: Target structure (left), posterior mean estimate (center) and
approximate standard deviations (right) of the target structure in the
simulated label fusion example. (The color scale in the right plot is
different from the other two for better definition.)

Most existing label fusion procedures result in a binary map corresponding to
inclusion/exclusion of voxels in a structure of interest. Thresholding the
posterior probability maps produced through our proposed approach also yields
such a binary map. A reasonable threshold is 0.5. Figure 4 displays the
thresholded posterior probability map. For comparison, the Figure displays
also the segmentations that result from simple majority voting, globally-
weighted majority voting (Artaechevarria et al., 2008) in which each atlas is
weighted by the inverse of its average squared intensity difference, and
locally-weighted majority voting (Artaechevarria et al., 2009) in which the
contribution of each atlas at each pixel is weighted by the inverse squared
intensity difference from the target. We can see the superior recovery of the
target under our approach versus different versions of majority voting. The
Figure also gives the Dice similarity coefficients quantifying similarity
between each segmentation and the target, where the target is dichotomized to
zero / non-zero pixels. Letting $d_{v}$ denote the indicator of inclusion of
pixel $v$ in the estimated image, the Dice coefficient is calculated as
$D:=2\sum_{v=1}^{V}d_{v}T_{v}/(2\sum_{v=1}^{V}d_{v}T_{v}+\sum_{v=1}^{V}d_{v}(1-T_{v})+\sum_{v=1}^{V}(1-d_{v})T_{v})$,
with values close to one indicating strong agreement between two images, and
values close to zero indicating strong disagreement. Thresholding the
posterior probability map has resulted in $175\%$ improvement in similarity
over even the most favorable version of majority voting.

Figure 4: Left four panels: Binary segmentation estimates and Dice coefficents
($D$) obtained via thresholded posterior probability map, simple majority
voting, globally weighted (GW) majority voting, and locally weighted (LW)
majority voting in the simulated label fusion example. Right panel: Marginal
posterior density of volume along with the competing volume estimates.

Often, a researcher is uninterested in a binary segmentation for its own sake,
but rather as a means to an end; e.g., volume estimation. With a binary
segmentation, the only way to estimate the volume is by summing the binary
indicators over the image. In contrast, our approach facilitates construction
of a distribution of plausible volumes. Let $T^{(k)}_{v}$ denote the
realization of $T_{v}$ on iteration $k$ of the MCMC output. Recognizing that
$N^{-1}\sum_{k=1}^{N}\sum_{v=1}^{V}T_{v}^{(k)}\approx\sum_{v=1}^{V}E(T_{v}\mid\boldsymbol{Y})=\sum_{v=1}^{V}E[P(T_{v}=1\mid\boldsymbol{\omega},\boldsymbol{Y})]$,
we estimate the structure volume on iteration $k$ with
$\hat{M}^{(k)}:=\sum_{v=1}^{V}\hat{P}(T_{v}=1\mid\boldsymbol{\omega}^{(k)},\boldsymbol{Y}),~{}k=1,\ldots,N$.
The right panel in Figure 4 displays the estimated marginal distribution of
volume, along with the true volume and the volume estimates obtained from the
three voting procedures. The true volume is well within the high probability
region of the distribution. We have no way of formally quantifying the
uncertainty associated with the voting estimates.

These simulation results demonstrate possible improvements over simple and
weighted majority voting through our proposed Bayesian label fusion model.
This is possible even when strongly corrupted atlases are used as inputs to
the label fusion, since covariate and prior information can protect against
otherwise poor segmentaitons. By allowing each atlas’ reliability to vary
throughout an image, we account for the fact that no atlas is uniformly more
reliable than another throughout the image.

##  4 Application to Hippocampus Segmentation

The Alzheimer’s Disease Neuroimaging Initiative (ADNI;
http://www.adni.loni.usc.edu) is an ongoing, multi-center, longitudinal study
with the goals of better understanding Alzheimer’s disease (AD) and developing
effective biomarkers. It is important to gain a better understanding of early
disease mechanisms in the brain to inspire the development of effective,
targeted therapies for the disease. The association between hippocampal
atrophy and AD is well-documented (Likeman et al., 2005; Leung et al., 2010)
and volumetric changes over time can be used to aid the diagnosis of early AD
and track progression (Dubois et al., 2007).

Using demographic information available from the ADNI, we created an age- and
sex-matched, case-control sample of six AD and six healthy control subjects.
Manual segmentations of the left and right hippocampus were obtained for each
subject from the Harmonized Protocol For Hippocampal Volumetry (Boccardi et
al., 2015). The corresponding T1 images were downloaded from the LONI IDA
website. For each image, we applied N4 inhomogeneity correction (Tustison et
al., 2010) and Multi-atlas Skull Stripping (Doshi et al., 2013). For each
target T1 image we consider, the remaining T1 atlases were non-linearly
registered to the target image using `SyN` deformable registration with mutual
information cost and Welch windowed sinc interpolation (Avants et al., 2011a)
in `R`. The concomitant transformations were applied to each manual
segmentation with nearest-neighbor interpolation to obtain the atlases for
each target. The T1 intensity values were normalized across images using white
stripe (Shinohara et al., 2014). A single axial slice was chosen for each
target image that approximately maximized the hippocampal volume for that
subject across slices.

In this Section, we are motivated by a researcher using known control
subjects’ brain images as the atlas set for segmenting both healthy brains and
those that have been diagnosed with AD. We first display results from using
the control images as atlases for segmenting a single brain image that has
been diagnosed with AD. We subsequently compare our proposed method to
standard label fusion methods over all six AD brain images. In addition, we
summarize performance on healthy brains by separately considering each control
subject as a target image to be segmented with with remaining control subjects
as atlases.

###  4.1 Segmenting the Hippocampus of a Single Diseased Brain

We take as our target the brain image of a 78 year old male diagnosed with AD.
The six controls are registered to the target brain, resulting in the atlas
segmentations displayed in Supplementary Figure 4. As a covariate, we take the
weighted sum of the SDL maps for each atlas (rescaled to the unit interval),
where each atlas is weighted voxelwise by its T1 intensity similarity with the
target T1 image. As the hippocampus is a structure consisting of gray matter,
our second covariate is a binary gray matter indicator for each voxel. The
tissue class segmentation is obtained via the ATROPOS algorithm (Avants et
al., 2011b). The SDL map and tissue class segmentation, along with the
manually-segmented target image, are displayed in Figure 5. The target image
in this case has resolution $113\times 101$ so that $V=11,413$ and $R=6$ in
(1). We include an interaction term to downweight the effect of the presence
of gray matter when a voxel is thought to be far away from the hippocampus. To
elicit the prior on the regression coefficients via CMP, we use pseudodata to
impose our prior knowledge that (i) gray matter and a small SDL has a very
high likelihood of truly belonging to the hippocampus, (ii) a large SDL and
not gray matter is very likely to be external to the hippocampus, (iii) gray
matter with a large SDL is very likely external to the hippocampus, and (iv) a
small SDL but segmented as not belonging to gray matter has a marginal
probability of being part of the hippocampus, due to possibly misclassified
tissue type. (See Supplmentary Material for additional details.) The priors on
the sensitivity and specificity processes are specified as they were in
Section 3. After initialization, we run 100,000 iterations of a single Monte
Carlo Markov chain using the algorithm described in Subsection 2.3, thinning
to every $50^{\text{th}}$ draw. The last 1,000 realizations are retained as an
approximate sample from the posterior distribution. We examine trace plots and
ergodic averages of the $\boldsymbol{\delta}$ coefficients as a convergence
diagnostic. Supplementary Figure 5 displays these plots along with the Geweke
statistics and lag 1 autocorrelations, all of which suggest convergence. The
entire MCMC routine is coded in R and executed on a Dell Precision T3620
desktop PC running Windows 10 with a Intel Xeon 4.1 GHz CPU and 64GB RAM. It
takes approximately 2.3 hours to run.

Figure 5: Manually segmented hippocampus (left), weighted and summed SDL map
of the atlases (middle) and tissue class segmentation of the target (right)
for the label fusion application. In the tissue class segmentation, black,
dark gray, light gray, and white correspond to background, cerebrospinal
fluid, gray matter, and white matter, respectively.

Figure 6 displays the approximate posterior probability map for hippocampus
inclusion, along with the point-wise standard deviations. The highest
probabilities of inclusion are near the centers of the two regions. There is
considerable uncertainty about the edges. If a researcher were only interested
in the hippocampal volume for the patient, these probabilities could be used
directly to estimate the volume, as discussed below.

Figure 6: Posterior inclusion probability map and pointwise standard
deviations for segmenting the hippocampus of an AD patient.

One can threshold the posterior probabilities to obtain a binary inclusion
map. After thresholding, we compare the resulting segmentation to those
obtained by simple majority voting, globally-weighted majority voting,
locally-weighted majority voting, and joint label fusion (JLF; Wang et al.,
2013). The global weighting is inversely proportional to each atlases average
T1 intensity difference from the target. Local weighting is done similarly
using voxel-specific intensity differences. Our implementation of JLF takes as
input the same two-dimensional slices as the other procedures for a direct
comparison.

Figure 7 displays the manual segmentation along with that which is obtained by
thresholding the posterior inclusion probabilities at 0.5. Bayesian label
fusion attains a Dice coefficient of $D=0.773$. The next closest segmentation
is produced by JLF ($D=0.693$), followed by locally-weighted majority voting
($D=0.647$), simple majority voting ($D=0.632$) and globally-weighted majority
voting ($D=0.631$). Since only healthy brains are used as atlases, there is a
clear tendency for the established methods to over-segment the hippocampus. In
contrast, by explicitly incorporating the estimated gray matter pattern as a
predictor, our model is able to more closely follow the visible contours in
the image. The result is considerably stronger agreement between our approach
and the manual segmentation.

Figure 7: Manual (top left) and automatic segmentations produced by the
proposed Bayesian label fusion (BLF, top middle) model, joint label fusion
(JLF, top right), locally-weighted (LW, bottom left) majority voting,
globally-weighted (GW, bottom middle) majority voting, and simple majority
voting (bottom right). Displayed above each automated segmentation is the Dice
similarity coefficient ($D$) and absolute volume difference ($AVD$) between
the automatic and manual segmentations.

In practice, an anatomical structure is segmented to obtain important
information such as its volume or average image intensity within the
structure. In our case, segmenting the hippocampus is a step toward estimating
its volume. If we only obtain a binary map, then the only way to estimate the
volume is by summing the indicators. Doing so ignores many sources of
uncertainty, including image pre-processing, registration error, biological
variation, and rater variability. Additional uncertainty arises from partial
volume effects in which a voxel may not be entirely inside or outside a
structure of interest, but partly in both. By providing probabilities of
inclusion, our results can be interpreted as partly accounting for partial
volume effects (Ashburner and Friston, 2005). Monte Carlo sampling also
facilitates estimation of a distribution of plausible volumes through
$\hat{M}^{(k)}$, defined in Section 3. Figure 8 displays the marginal volume
density for the diseased brain of interest. Vertical lines indicate the
manually-segmented volume and the volume estimates from the competing
procedures. Since JLF and the voting procedures produce only point estimates,
it is difficult to meaningfully associate uncertainty with these estimates.
There is also uncertainty associated with the manual segmentation, since a
radiologist is unlikely to produce identical segmentations of the same brain
twice. Nevertheless, the distributional volume estimate produced by our
approach agrees with the manual segmentation more so than the other
procedures. We quantify the difference in volume estimates with the absolute
volume difference between the manual volume estimate, $\hat{M}_{m}$, and an
automatically segmented volume estimate, $\hat{M}_{a}$, (Carass et al., 2017),
defined as
$AVD(\hat{M}_{a},\hat{M}_{m})=(\max\\{\hat{M}_{a},\hat{M}_{m}\\}-\min\\{\hat{M}_{a},\hat{M}_{m}\\})/\hat{M}_{m}$,
where smaller values indicate stronger agreement. Figure 7 also displays the
AVDs for each automatic segmentation. Our proposed procedure produces a much
smaller difference ($AVD=0.160$) than JLF ($AVD=0.639$), locally-weighted
majority voting ($AVD=0.528$), globally-weighted majority voting ($AVD=0.631$)
or simple majority voting ($AVD=0.528$).

Figure 8: Estimated marginal posterior distribution of intraslice hippocampus
volume produced by the proposed Baysian label fusion model. Also plotted are
the manual, joint label fusion (JLF), locally-weighted (LW), globally-weighted
(GW), and simple majority voting volume estimates. (In this example, LW and
simple majority voting produce identical volume estimates.)

Our proposed model also has the ability to estimate the spatially-varying
reliability parameters for each atlas. Figure 9 displays such reliability maps
for two selected atlases. We see the smooth decay of sensitivity between
regions of high agreement and low agreement with the manual segmentations.
Since the sensitivity is conditional on the voxel truly being part of the
structure of interest ($T_{v}=1$), the model is only capable of estimating
sensitivity in areas where it estimates a high probability of the voxel truly
being part of the structure. This is evident in the sensitivity maps. As we
move away from the structure estimate, there is no information with which to
update the spatial fields. Thus the fields return to the prior mean.
Similarly, the specificity values are estimated to be very high away from the
targeted structure, where all the atlases agree on the voxels being excluded.
The specificity fields are determined by their prior mean in areas where
$P(T_{v}=0\mid\boldsymbol{y})$ is quite low.

Figure 9: Estimated reliability fields for two selected atlases used in
segmenting the hippocampus of one diseased brain from the ADNI data. Also
displayed are the manual segmentation and the actual atlas segmentation for
reference.

###  4.2 Aggregated Results

To study how our proposed approach performs more generally, we present the
results of using the control brain images as atlases for each of the six AD
patients in our dataset. The top row of Figure 10 summarizes the Dice
coefficients, volume estimates, and absolute volume differences between the
automatic and manual segmentations. Supplementary Figure 6 displays manual
segmentations of two additional AD patients along with the posterior
probability and pointwise standard deviation maps from our approach. In terms
of image similarity, our approach and JLF separate themselves from the three
voting approaches. The top middle panel shows that the volume estimates from
JLF are consistently greater than the manual volume estimates, whereas the
posterior mean volume estimates produced by our procedure tend to be closer.
The result is that the absolute volume differences from our method tend to be
much smaller than those from the other procedures. In addition, out of any of
the volume estimates (including manual segmentation), ours is the only one
that meaningfully produces measures of associated uncertainty, depicted in
this case as 99% posterior credible intervals. Bayesian posterior
distributions reflect uncertainty about parameters in the presence of the
observed data balanced against a priori knowledge under an assumed statistical
model; a model which is at best only an approximation to reality. Thus,
posterior credible intervals do not necessarily exhibit frequentist operating
characteristics, nor should they be expected to. Even a trained expert will
produce slightly different manual segmentations of the brain on two different
occasions, so there is uncertainty associated with each manual segmentation
itself. This is not quantified here. Thus it is impossible to assess any
‘significant’ difference between volume estimates from our proposed approach
and a manually estimated volume.

Figure 10: Summary of Dice similarity coefficients, volume estimates, and
absolute volume differences for the proposed Bayesian label fusion (BLF),
joint label fusion (JLF), locally-weighted majority voting (LWMV), globally-
weighted majority voting (GWMV), and simple majority voting (MV) when using
control cases as atlases for the AD patients. Both the posterior mean (BLF_PM)
and thresholded probabilities (BLF_TH) are used as volume estimates under the
BLF approach. The bars around the BLF estimates in the middle panel depict
approximate 99% credible intervals.

Supplementary Figure 7 displays the results from an experiment in which each
of the six control images is held out and segmented using the remaining five
controls as atlases. There is much less distinction between the competing
approaches, though JLF tends to produce the highest Dice coefficients. Our
approach produces Dice coefficients that are much more variable compared to
the AD scenario, though they are still mostly above 0.8. In terms of volume
differences, all of the approaches are comparable. This is not surprising
since using healthy brains to segment healthy brains is an ideal case when
compared to the more challenging task of segmenting diseased brains.

Our assumed regression model is the same for both the diseased and control
scenarios. This suggests that relying on auxiliary tissue class information
can be quite helpful when segmenting diseased brains, whereas such information
may not be as important in the healthy case. Our proposed approach with
different image intensity-derived covariates and/or different prior
elicitations could possibly produce stronger results for the control cases. We
defer this issue to future investigation.

Application of our proposed model to label fusion in the ADNI data
demonstrates the feasibility of our approach as well as additional useful
information that would be unavailable otherwise. Posterior probability maps
can be used to summarize where the hippocampus is likely to be. A simple
thresholding rule yields image similarities that outperform even state of the
art methods when using healthy brains as atlases for diseased cases.
Thresholding is unnecessary when the goal is to estimate the volume of the
hippocampus, as shown by the availability of a posterior volume distribution.

##  5 Discussion

Segmenting brain images is a challenging task that becomes quite resource-
intensive when done manually. Nevertheless, it is necessary for
neuroscientists and neurologists to accurately measure traits that can be
related to clinical endpoints. Thus there is a need for the continued
development of automated procedures for segmentation. The utility of
segmentation is not limited to structural brain images. Functional MRI also
frequently uses ROIs (e.g., Bowman, 2007). While ROI definitions in functional
MRI sometimes differ from anatomical definitions, the problem of delineating
those regions in a brain remains the same.

We have taken a hierarchical Bayesian approach to label fusion for anatomical
segmentation. We proposed a spatial binary regression model that not only
determines the likely location of the target structure, but estimates the
reliability of each individual atlas via spatially-varying sensitivity and
specificity fields. These reliability parameters are built in to the Bayesian
model so that each atlas’ ‘vote’ is automatically weighted when estimating the
inclusion probability at each voxel. We were able to fully account for all
sources of uncertainty by accessing the posterior distribution via MCMC. With
posterior probability maps and their associated distributions, a researcher
can threshold to obtain a binary segementation. Our proposed approach also
facilitates estimating a distribution of plausible volumes.

We considered in this work segmentation of two-dimensional slices from
multiple brain images. We believe the proposed approach can extend to three-
dimensional, multi-slice segmentation of hippocampus with little difficulty.
In fact, the new ADNI 3 studies include hippocampal scans consisting of six to
seven slices, and chromatic sampling has been shown to scale at a much slower
rate than block sampling (Brown et al., 2019). In some cases, many more slices
are needed, either for larger binary segmentation problems or whole-brain
parcellation. Toward this end, there exists in the literature a variety of
techniques for MCMC in high dimensional spaces, including when the model
involves GMRFs (e.g., Sidén et al., 2017). Additional options are available if
one is willing to move away from GMRFs. For instance, Bezener et al. (2018)
use pre-defined regions to aggregate voxels and reduce the dimension of the
spatial field underlying fMRI data while maintaining spatial dependence. Zhao
et al. (2018) mitigate the computational burden with a multi-resolution MCMC
approach that successively refines resolutions in interesting areas of a brain
image. Heaton et al. (2018) review many large scale Gaussian process
techniques that have been proposed recently. Also available are posterior
approximation methods such as variational Bayes, though such approximations
can yield inferior results to those produce by MCMC (Teng et al., 2016; Sidén
et al., 2017). Concerning the extension from binary segmentation to whole-
brain parcellation, an obvious modification to our proposed approach is to
replace the binomial response with multinomial. This raises new challenges,
not the least of which is computation. Given the aformentioned MCMC advances,
though, we are optimistic that such a path is within reach.

##  6 Acknowledgements

This material is based upon work partially supported by the National Science
Foundation (NSF) under Grant DMS-1127914 to the Statistical and Applied
Mathematical Sciences Institute. DAB is partially supported by NSF Grants
CMMI-1563435, EEC-1744497, and OIA-1826715. CSM is partially supported by
National Institutes of Health Grant AI121351 and NSF grant OIA-1826715. RTS is
partially supported by NIH Grants R01NS085211, R21NS093349, R01MH112847,
R01NS060910, and R01EB017255. The content is solely the responsibility of the
authors and does not necessarily represent the official views of the funding
agencies. The authors thank the Editor, an Associate Editor, anonymous
referees, and the clinical neuroimaging research group in the Department of
Biostatistics at Johns Hopkins University, particularly John Muschelli for his
help with the image registrations.

## References

  * Akhondi-Asl et al. (2014) Akhondi-Asl, A., Hoyte, L., Lockhart, M., and Warfield, S. (2014).  A logarithmic opinion pool based STAPLE algorithm for the fusion of segmentation with associated reliability weights.  IEEE Transactions on Medical Imaging, 33(10):1997–2009.
  * Akkus et al. (2017) Akkus, Z., Galimzianova, A., Hoogi, A., Rubin, D. L., and Erickson, B. L. (2017).  Deep learning for brain MRI segmentation: State of the art and future directions.  Journal of Digitial Imaging, 30:449–459.
  * Albert and Chib (1993) Albert, J. H. and Chib, S. (1993).  Bayesian analysis of binary and polychotomous data.  Journal of the American Statistical Association, 88(422):669–679.
  * Artaechevarria et al. (2008) Artaechevarria, X., Muñoz Barrutia, A., and Ortiz-de Solorzano, C. (2008).  Efficient classifier generation and weighted voting for atlas-based segmentation: Two small steps faster and closer to the combination oracle.  In Medical Imaging, page 69141W. International Society for Optics and Photonics.
  * Artaechevarria et al. (2009) Artaechevarria, X., Muñoz Barrutia, A., and Ortiz-de Solorzano, C. (2009).  Combination strategies in multi-atlas image segmentation: Application to brain MR data.  IEEE Transactions on Medical Imaging, 28(8):1266–1277.
  * Ashburner and Friston (2005) Ashburner, J. and Friston, K. J. (2005).  Unified segmentation.  NeuroImage, 26(3):839–851.
  * Avants et al. (2011a) Avants, B. B., Tustison, N. J., Song, G., Cook, P. A., Klein, A., and Gee, J. C. (2011a).  A reproducible evaluation of ANTs similarity metric performance in brain image registration.  NeuroImage, 54(3):2033–2044.
  * Avants et al. (2011b) Avants, B. B., Tustison, N. J., Wu, J., Cook, P. A., and Gee, J. C. (2011b).  An open source multivariate framework for n-tissue segmentation with evaluation on public data.  Neuroinformatics, 9(4):381–400.
  * Banerjee et al. (2015) Banerjee, S., Carlin, B. P., and Gelfand, A. E. (2015).  Hierarchical Modeling and Analysis for Spatial Data.  Chapman & Hall/CRC, Boca Raton, 2nd edition.
  * Bedrick et al. (1996) Bedrick, E. J., Christensen, R., and Johnson, W. (1996).  A new perspective on priors for generalized linear models.  Journal of the American Statistical Association, 91(436):1450–1460.
  * Bernardinelli et al. (1995) Bernardinelli, L., Clayton, D. G., and Montomoli, C. (1995).  Bayesian estimates of disease maps: How important are priors?  Statistics in Medicine, 14:2411–2431.
  * Besag (1974) Besag, J. (1974).  Spatial interaction and the statistical analysis of lattice systems (with discussion).  Journal of the Royal Statistical Society, Series B, 36:192–236.
  * Besag and Kooperberg (1995) Besag, J. and Kooperberg, C. (1995).  On conditional and intrinsic autoregressions.  Biometrika, 82(4):733–746.
  * Bezener et al. (2018) Bezener, M., Hughes, J., and Jones, G. (2018).  Bayesian spatiotemporal modeling using hierarchical spatial priors, with applications to functional magnetic resonance imaging (with discussion).  Bayesian Analysis, 13(4):1261–1313.
  * Boccardi et al. (2015) Boccardi, M., Bocchetta, M., Morency, F. C., Collins, D. L., Nishikawa, M., Ganzola, R., Grothe, M. J., Wolf, D., Redolfi, A., Pievani, M., et al. (2015).  Training labels for hippocampal segmentation based on the EADC-ADNI harmonized hippocampal protocol.  Alzheimer’s and Dementia, 11(2):175–183.
  * Bowman (2007) Bowman, F. D. (2007).  Spatiotemporal models for region of interest analyses of functional neuroimaging data.  Journal of the American Statistical Association, 102(478):442–453.
  * Brown et al. (2019) Brown, D. A., McMahan, C. S., and Self, S. W. (2019).  Sampling strategies for fast updating of Gaussian Markov random fields.  ArXiv:1702.05518.
  * Carass et al. (2017) Carass, A., Roy, S., Jog, A., Cuzzocreo, J. L., Magrath, E., Gherman, A., Button, J., Nguyen, J., Prados, F., Sudre, C. H., et al. (2017).  Longitudinal multiple sclerosis lesion segmentation: Resource and challenge.  NeuroImage, 148:77–102.
  * Doshi et al. (2013) Doshi, J., Erus, G., Ou, Y., Gaonkar, B., and Davatzikos, C. (2013).  Multi-atlas skull-stripping.  Academic Radiology, 20(12):1566–1576.
  * Doshi et al. (2016) Doshi, J., Erus, G., Ou, Y., Resnick, S. M., Gur, R. C., Gur, R. E., Satterthwaite, T. D., Furth, S., and Davatzikos, C. (2016).  MUSE: MUlti-atlas Segmentation utilizing Ensembles of registration algorithms and parameteres, and locally optimal atlas selection.  NeuroImage, 127:186–195.
  * Dubois et al. (2007) Dubois, B., Feldman, H. H., Jacova, C., DeKosky, S. T., Barberger-Gateau, P., Cummings, J., Delacourte, A., Galasko, D., Gauthier, S., Jicha, G., et al. (2007).  Research criteria for the diagnosis of Alzheimer’s disease: Revising the NINCDS–ADRDA criteria.  The Lancet Neurology, 6(8):734–746.
  * Eberly and Carlin (2000) Eberly, L. E. and Carlin, B. P. (2000).  Identifiability and convergence issues for Markov chain Monte Carlo fitting of spatial models.  Statistics in Medicine, 19:2279–2294.
  * Fox and Nicholls (2001) Fox, C. and Nicholls, G. K. (2001).  Exact MAP states and expectations from perfect sampling: Greig, Porteous and Seheult revisited.  In Mohammed-Djafari, A., editor, AIP Conference Proceedings, volume 568, pages 252–263. AIP.
  * Gamerman (1997) Gamerman, D. (1997).  Sampling from the posterior distribution in generalized linear mixed models.  Statistics and Computing, 7:57–68.
  * Gelfand and Smith (1990) Gelfand, A. E. and Smith, A. F. M. (1990).  Sampling-based approaches to calculating marginal densities.  Journal of the American Statistical Association, 85(410):398–409.
  * Geman and Geman (1984) Geman, S. and Geman, D. (1984).  Stochastic relaxation, Gibbs distributions and the Bayesian restoration of images.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 6:721–741.
  * Geweke (1992) Geweke, J. (1992).  Evaluating the accuracy of sampling-based approaches to calculation of posterior moments.  In Proceedings of the Fourth Valencia International Conference on Bayesian Statistics, pages 169–193.
  * Gonzalez et al. (2011) Gonzalez, J. E., Low, Y., Gretton, A., and Guestrin, C. (2011).  Parallel Gibbs sampling: From colored fields to thin junction trees.  In Journal of Machine Learning Research: Proceedings of the 14th International Conference on Artificial Intelligence and Statistics (AISTATS), pages 324–332.
  * Greve and Fischl (2009) Greve, D. N. and Fischl, B. (2009).  Accurate and robust brain image alignment using boundary-based registration.  NeuroImage, 48(1):63–72.
  * Hanson et al. (2012) Hanson, J. L., Suh, J. W., Nacewicz, B. M., Sutterer, M. J., Cayo, A. A., Stodola, D. E., Burghy, C. A., Wang, H., Avants, B. B., Yushkevich, P. A., et al. (2012).  Robust automated amygdala segmentation via multi-atlas diffeomorphic registration.  Frontiers in Neuroscience, 6:1–8.
  * Hastings (1970) Hastings, W. (1970).  Monte Carlo sampling methods using Markov chains and their application.  Biometrika, 57:97–109.
  * Heaton et al. (2018) Heaton, M. J., Datta, A., Finley, A. O., Furrer, R., Guinness, J., Guhaniyogi, R., Gerber, R., Gramacy, R. B., Hammerling, D., Katzfuss, M., Lindgren, F., Nychka, D. W., Sun, F., and Zammit-Mangion, A. (2018).  A case study competition among methods for analyzing large spatial data.  Journal of Agricultural, Biological, and Environmental Statistics.  To appear.
  * Heckemann et al. (2006) Heckemann, R. A., Hajnal, J. V., Aljabar, P., Rueckert, D., and Hammers, A. (2006).  Automatic anatomical brain MRI segmentation combining label propagation and decision fusion.  NeuroImage, 33(1):115–126.
  * Huo et al. (2016) Huo, Y., Plassard, A. J., Carass, A., Resnick, S. M., Pham, D. L., Prince, J. L., and Landman, B. A. (2016).  Consistent cortical reconstruction and multi-atlas brain segmentation.  NeuroImage, 138:197–210.
  * Iglesias and Sabuncu (2015) Iglesias, J. E. and Sabuncu, M. R. (2015).  Multi-atlas segmentation of biomedical images: A survey.  Medical Image Analysis, 24(1):205–219.
  * Iglesias et al. (2012) Iglesias, J. E., Sabuncu, M. R., and Leemput, K. V. (2012).  A generative model for probabilistic label fusion of multimodal data.  In Multimodal Brain Image Analysis, pages 115–133.
  * Iglesias et al. (2013) Iglesias, J. E., Sabuncu, M. R., and Leemput, K. V. (2013).  Improved inference in Bayesian segmentation using Monte Carlo sampling: Application to hippocampal subfield volumetry.  Medical Image Analysis, 17:766–778.
  * Isgum et al. (2009) Isgum, I., Staring, M., Rutten, A., Prokop, M., Viergever, M. A., and van Ginneken, B. (2009).  Multi-atlas-based segmentation with local decision fusion: Application to cardiac and aortic segmentation in CT scans.  IEEE Transactions in Medical Imaging, 28(7):1000–1010.
  * Jenkinson et al. (2002) Jenkinson, M., Bannister, P., Brady, J. M., and Smith, S. M. (2002).  Improved optimisation for the robust and accurate linear registration and motion correction of brain images.  NeuroImage, 17(2):825–841.
  * Jenkinson et al. (2012) Jenkinson, M., Beckmann, C. F., Behrens, T. E., Woolrich, M. W., and Smith, S. M. (2012).  FSL.  NeuroImage, 62:782–790.
  * Jenkinson and Smith (2001) Jenkinson, M. and Smith, S. M. (2001).  A global optimisation method for robust affine registration of brain images.  Medical Image Analysis, 5(2):143–156.
  * Klein et al. (2005) Klein, A., Mensh, B., Ghosh, S., Tourville, J., and Hirsch, J. (2005).  Mindboggle: Automated brain labeling with multiple atlases.  BMC Medical Imaging, 5(1):7.
  * Langerak et al. (2010) Langerak, T. R., van der Heide, U. A., Kotte, A. N. T. J., Viergever, M. A., van Vulpen, M., and Pluim, J. P. W. (2010).  Label fusion in atlas-based segmentation using a selective and iterative method for performance level estimation (SIMPLE).  IEEE Transactions on Medical Imaging, 29(12):2000–2008.
  * Leung et al. (2010) Leung, K. K., Barnes, J., Ridgway, G. R., Bartlett, J. W., Clarkson, M. J., Macdonald, K., Schuff, N., Fox, N. C., Ourselin, S., et al. (2010).  Automated cross-sectional and longitudinal hippocampal volume measurement in mild cognitive impairment and Alzheimer’s disease.  NeuroImage, 51(4):1345–1359.
  * Likeman et al. (2005) Likeman, M., Anderson, V. M., Stevens, J. M., Waldman, A. D., Godbolt, A. K., Frost, C., Rossor, M. N., and Fox, N. C. (2005).  Visual assessment of atrophy on magnetic resonance imaging in the diagnosis of pathologically confirmed young-onset dementias.  Archives of Neurology, 62(9):1410–1415.
  * Metropolis et al. (1953) Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. H., and Teller, E. (1953).  Equation of state calculations by fast computing machines.  Journal of Chemical Physics, 21:1087–1091.
  * Milletari et al. (2017) Milletari, F., Ahmadi, S.-A., Kroll, C., Plate, A., Rozanski, V., Maiostre, J., Levin, J., Dietrich, O., Ertl-Wagner, B., Bötzel, K., and Navab, N. (2017).  Hough-CNN: Deep learning for segmentation of deep brain regions in MRI and ultrasound.  Computer Vision and Image Understanding, 164:92–102.
  * Müller (1991) Müller, P. (1991).  A generic approach to posterior integration and Gibbs sampling.  Technical Report, Purdue University.
  * Ou et al. (2011) Ou, Y., Sotiras, A., Paragios, N., and Davatzikos, C. (2011).  DRAMMS: Deformable registration vis attribute matching and mutual-saliency weighting.  Medical Image Analysis, 15(4):622–639.
  * Park et al. (2014) Park, M. T. M., Pipitone, J., Bear, J. H., Winterburn, J. L., Shah, Y., Chavez, S., Schira, M. M., Lobaugh, N. J., Lerch, J. P., Voineskos, A. N., et al. (2014).  Derivation of high-resolution MRI atlases of the human cerebellum at 3T and segmentation using multiple automatically generated templates.  NeuroImage, 95:217–231.
  * Pham et al. (2000) Pham, D. L., C., X., and Prince, J. L. (2000).  Current methods in medical image segmentation 1.  Annual Review of Biomedical Engineering, 2(1):315–337.
  * R Development Core Team (2008) R Development Core Team (2008).  R: A Language and Environment for Statistical Computing.  R Foundation for Statistical Computing, Vienna, Austria.  ISBN 3-900051-07-0.
  * Raamana et al. (2014) Raamana, P. R., Rosen, H., Miller, B., Weiner, M. W., Wang, L., and Beg, M. F. (2014).  Three-class differential diagnosis among Alzheimer’s disease, frontotemporal dementia, and controls.  Frontiers in Neurology, 5.
  * Ramus and Malandain (2010) Ramus, L. and Malandain, G. (2010).  Assessing selection methods in the context of multi-atlas based segmentation.  In Proceedings of the IEEE International Symposium on Biomedical Imaging: From Nano to Macro, pages 1321–1324.
  * Rohlfing et al. (2004) Rohlfing, T., Brandt, R., Menzel, R., and Maurer Jr., C. R. (2004).  Evaluation of atlas selection strategies for atlas-based segmentation with application to confocal microscopy of bee brains.  NeuroImage, 21(4):1428–1442.
  * Rue and Held (2005) Rue, H. and Held, L. (2005).  Gaussian Markov Random Fields.  Chapman & Hall/CRC, Boca Raton.
  * Sabuncu et al. (2010) Sabuncu, M. R., Yeo, B. T., van Leemput, K., Fischl, B., and Golland, P. (2010).  A generative model for image segmentation based on label fusion.  IEEE Transactions on Medical Imaging, 29(10):1714–1729.
  * Shinohara et al. (2014) Shinohara, R. T., Sweeney, E. M., Goldsmith, J., Shiee, N., Mateen, F. J., Calabresi, P. A., Jarso, S., Pham, D. L., Reich, D. S., and Crainiceanu, C. M. (2014).  Statistical normalization techniques for magnetic resonance imaging.  NeuroImage: Clinical, 6:9–19.
  * Sidén et al. (2017) Sidén, P., Eklund, A., Bolin, D., and Villani, M. (2017).  Fast Bayesian whole-brain fMRI analysis with spatial 3D priors.  NeuroImage, 146:211–225.
  * Teng et al. (2016) Teng, M., Johnson, T., and Nathoo, F. s. (2016).  A comparison of variational Bayes and Hamiltonian Monte Carlo for Bayesian fMRI time series analysis with spatial priors.  ArXiv:1609.02123.
  * Tustison and Avants (2013) Tustison, N. J. and Avants, B. B. (2013).  Explicit B-spline regularization in diffeomorphic image registrations.  Frontiers in Neuroinformatics, 7.
  * Tustison et al. (2010) Tustison, N. J., Avants, B. B., Cook, P. A., Zheng, Y., Egan, A., Yushkevich, P. A., and Gee, J. C. (2010).  N4ITK: Improved N3 bias correction.  IEEE Transactions on Medical Imaging, 29(6):1310–1320.
  * Wachinger et al. (2018) Wachinger, C., Reuter, M., and Klein, T. (2018).  DeepNAT: Deep convolutional neural network for segmenting neuroanatomy.  NeuroImage, 170:434–445.
  * Wang et al. (2013) Wang, H., Suh, J. W., Sandhitsu, R. D., Pluta, J. B., Craige, C., and Yushkevich, P. A. (2013).  Multi-atlas segmentation with joint label fusion.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 35(3):611–623.
  * Wang et al. (2012) Wang, Y., Jia, H., Yap, P.-T., Cheng, B., Wee, C.-Y., Guo, L., and Shen, D. (2012).  Groupwise segmentation improves neuroimaging classification accuracy.  In Multimodal Brain Image Analysis, pages 185–193. Springer.
  * Warfield et al. (2004) Warfield, S. K., Zou, K. H., and Wells, W. M. (2004).  Simultaneous truth and performance level estimation (STAPLE): An algorithm for the validation of image segmentation.  IEEE Transactions on Medical Imaging, 23(7):903–921.
  * Yushkevich et al. (2006) Yushkevich, P. A., Piven, J., Hazlett, H. C., Smith, R. G., Ho, S., Gee, J. C., and Gerig, G. (2006).  User-guided 3D active contour segmentation of anatomical structures: Significantly improved efficiency and reliability.  NeuroImage, 31(3):1116–1128.
  * Zhao et al. (2018) Zhao, Y., Kang, J., and Long, Q. (2018).  Bayesian multiresolution variable selection for ultra-high dimensional neuroimaging data.  IEEE/ACM Transactions on Computational Biology and Bioinformatics, 15(2):537–550.
  * Zhong et al. (2010) Zhong, J., Phua, D. Y. L., and Qiu, A. (2010).  Quantitative evaluation of LDDMM, FreeSurfer, and CARET for cortical surface mapping.  NeuroImage, 52(1):131–141.