# Bayesian Functional Registration of fMRI Activation Maps

Guoqing Wang  Abhirup Datta  Martin A<EMAIL_ADDRESS>[
Department of Biostatistics, Johns Hopkins University,

###### Abstract

Functional magnetic resonance imaging (fMRI) has provided invaluable insight
into our understanding of human behavior. However, large inter-individual
differences in both brain anatomy and functional localization _after_
anatomical alignment remain a major limitation in conducting group analyses
and performing population level inference. This paper addresses this problem
by developing and validating a new computational technique for reducing
misalignment across individuals in functional brain systems by spatially
transforming each subjects functional data to a common reference map. Our
proposed Bayesian functional registration approach allows us to assess
differences in brain function across subjects and individual differences in
activation topology. It combines intensity-based and feature-based information
into an integrated framework, and allows inference to be performed on the
transformation via the posterior samples. We evaluate the method in a
simulation study and apply it to data from a study of thermal pain. We find
that the proposed approach provides increased sensitivity for group-level
inference.

functional magnetic resonance imaging,

group-level analysis,

registration,

Bayesian methods,

inter-individual differences,

pain,

###### keywords:

, and

##  1 Introduction

In recent years functional magnetic resonance imaging (fMRI) has provided
invaluable insight into understanding the neurophysiological underpinnings of
human behavior. Much of these gains have been driven by the combination of
increasingly sophisticated statistical and computational techniques and the
ability to capture brain data at finer spatial and temporal resolution.
Together, these advances are allowing researchers to study brain-behavior
correlations and develop population-level models of the functional brain
representations underlying behavior, performance, clinical status and
prognosis. However, a fundamental obstacle to using emerging analytic
techniques effectively in these settings is individual variation in functional
brain anatomy. That is, the fact that the precise locations of brain
activation are not perfectly aligned across subjects. Addressing this obstacle
will greatly extend the capability to perform group-level hypothesis testing,
use machine learning to develop population-level brain models that predict
behavior, and take advantage of the full spatial resolution of fMRI.

In a task-based fMRI study, each subject is typically administered one or more
stimuli and observed at hundreds of time points. At each time point, the
subject’s blood oxygenation level dependent (BOLD) response is recorded at
roughly 100,000 spatial locations (voxels), giving rise to multivariate time
series data. In order to properly perform statistical analysis over multiple
subjects, it is necessary for each voxel to lie within the same brain
structure for each subject. This does not occur naturally, as all subjects
have brains of different sizes and shapes. Therefore, to achieve this goal all
subjects’ brains are normalized to a stereotaxic space prior to analysis.

Typically, normalization is performed using a high-resolution anatomical scan
that is acquired in the same session and spatially aligned with the fMRI data.
Current standard practice is to use nonlinear transformations to warp
individual subjects’ anatomical scan to an average, anatomically based 3D
reference space (e.g., the “Montreal Neurologic Institute” (MNI) space). These
transformations are then applied to the fMRI data to ensure they lie in MNI
space, and data can be compared across subjects (Lindquist et al.,, 2008;
Ombao et al.,, 2016).

While this procedure efficiently makes the brains of all subjects overlap with
regards to gross anatomical landmarks, it does not implicitly take into
consideration residual differences in brain anatomy, and the location and
distribution of functional regions relative to these landmarks. For example,
primary visual cortex can vary in size by as much as 2-fold across different
subjects’ brains (Rademacher et al.,, 1993) and in location relative to other
anatomical landmarks (Amunts et al.,, 2000), as can sulcal locations (Thompson
et al.,, 1996). In some cases, inter-individual differences can be accounted
for by normalization, but in other cases they cannot. For example, the
cingulate sulcus is structurally dimorphic (Vogt et al.,, 1995) —a substantial
minority of individuals possesses a double cingulate sulcus—precluding
anatomical alignment using current methods. In addition, considerable
variability in functional localization persists even after anatomical
alignment (Duncan et al.,, 2009). For example, the location of visual motion-
sensitive area MT can vary by more than $2\,cm$. The location and extent of
lateral occipital cortex (Iordan et al.,, 2016) and fusiform face area
(Allison et al.,, 1999; McCarthy et al.,, 1999) are also highly variable
across individuals and located in different places relative to anatomical
sulcal landmarks.

Large inter-individual differences in both brain anatomy and functional
localization after anatomical alignment are a major limitation in conducting
group analyses and population level inference. This is true because a large
number of statistical models are performed voxel-wise (i.e., subjects are
compared using a separate model at each voxel). This can lead to a situation
where voxels from different functional brain regions are compared to one
another in a model, significantly complicating interpretation. In addition,
voxels are often used as features in prediction models, and any variation in
the location of functional regions across subjects would lead to feature
misalignment. Finally, it degrades the quality of the data. Though individual-
subject, whole-brain fMRI data can now be collected at a resolution of
$2\,mm^{3}$ voxels or smaller, once functional maps are combined across
individuals to construct a population-level model, the effective spatial
resolution is dramatically lower.

Recent approaches have sought to circumvent this problem by mapping individual
brains into a functional population-level reference space rather than an
anatomically based brain space. The first such model is the ‘hyperalignment’
procedure (Haxby et al.,, 2011). Here brain activity patterns corresponding to
stimuli and other cognitive events are represented as vectors in a neural
representational space spanned by the voxels in a local neighborhood.
Hyperalignment rotates each participant’s local voxel-wise activity patterns
through multivariate voxel space using a Procrustes transformation to align
the representational geometry across subjects. Mathematically, this is similar
to Canonical Correlation Analysis (CCA). A number of refinements of
hyperalignment have recently appeared, including: kernel hyperalignment
(Lorbert and Ramadge,, 2012), which performs nonlinear hyperalignment in an
embedding space; regularized hyperalignment (Xu et al.,, 2012) which uses a
ridge CCA formulation; the two-phase joint SVD-Hyperalignment algorithm (Chen
et al.,, 2014), where singular value decomposition (SVD) provides a lower
dimensional feature space where hyperalignment is applied; the shared response
model (Chen et al.,, 2014) which casts the model in a probabilistic framework;
and searchlight hyperalignment (Guntupalli et al.,, 2016) which allows for
whole-brain coverage. Other methods have been developed that align subjects
based on inter-subject correlations in time-series data during movie viewing.
One method is functional time series alignment (Sabuncu et al.,, 2010), which
matches voxels across subjects using a 2D ‘rubber sheet’ warping and
maximizing inter-subject correlations across voxels (Hasson et al.,, 2004).
Another is functional connectivity alignment (Conroy et al.,, 2013), which
matches voxels to a functional reference by minimizing the Frobenius norm of
the difference between a subject’s connectivity matrix and a reference matrix
with a shape-preserving penalty function. Both methods use warping and
penalization strategies specific to the cortex, so it is unclear how easily
they can be adapted to include subcortical regions.

This body of work shows great promise for allowing participants to vary in
their functional activation patterns (Haxby et al.,, 2011), but suffers from
several shortcomings. First, hyperalignment requires subjects to watch a long
film (up to 2 hours) to align subjects into a common representational space,
and functional connectivity-based methods require substantial resting state
data. Second, the choice of reference data influences which types of
functional patterns can be appropriately aligned. For example, movie reference
data works well for audio-visual representations but not prefrontal and limbic
networks. Functional connectivity-based measures work better for limbic cortex
but are expected to perform worse on object and semantic representations.
Neither type of method has been tested on subcortical representations, or on
clinically relevant functions such as pain and emotion. Third, these methods
do not include an explicit spatial model able to make inferences on the
location or extent of activation, though a few recent studies have taken steps
in this direction (Nenning et al.,, 2017).

In this work, we focus on performing local registration of brain activation
using subject-specific functional activation maps. These maps are the output
of a series of voxel-wise models that quantify the activation in the BOLD
signal that arises due to the stimulus of interest. They correspond to spatial
maps of the regression coefficient corresponding to the predictor modeling the
stimulus; for more detail see Section 2. As these are coefficient estimates
from a normal-error linear model that exhibit significant spatial correlation
they are well-described using a Gaussian process. It is important to note that
these models are fit after an initial anatomical registration has already been
performed. Thus, all subject-specific activation maps should be in a common
anatomical space, and all voxel-specific values directly comparable across
subjects. However, as discussed above, there remains substantial individual
variation in functional brain anatomy even after anatomical registration that
causes this assumption to be invalid. To circumvent this issue, we seek to
perform a secondary functional registration of the data into a common
functional space. This can be accomplished by warping (i.e., spatially
transforming) a set of floating maps (i.e., subject-specific activation maps)
to some predefined reference map (i.e., a particular subject-specific
activation map or a group-averaged map).

In service of this goal, we propose a generalized Bayes approach that is
capable of cropping the subject-specific activation maps to match certain pre-
defined regions of interest (ROIs). We adapt a flexible general loss function-
based pseudo-likelihood approach to measure and correct for relative
misalignment between the floating and reference maps. The pseudo-likelihood is
a representation of the similarity measure between the reference map and the
warped floating map. The measure can be the sum of squared difference (SSD) or
any other context-dependent choice. The parameters characterizing the
misalignment are assigned prior distributions estimated by matching the
features of the activation maps, which are commonly defined by local peaks.
Furthermore, our work simultaneously models an intensity correction term which
compensates for differences in the level of activation between the reference
and floating maps. As the images are defined on a regular lattice, the
discreteness leads to the need to utilize Gaussian process based Kriging
interpolation techniques (Banerjee et al.,, 2014; Cressie and Wikle,, 2015)
while realigning images. Finally, inference on the model parameters are based
on the draws from the posterior distribution. We use the Markov Chain Monte
Carlo (MCMC) algorithm to simulate from the posterior distributions.

This paper is organized as follows. Section 2 describes the data, which is
from an fMRI study of thermal pain. In Section 3 we introduce our generalized
Bayes framework for performing functional image registration. In Section 4 we
describe a simulation study used to evaluate the performance of the method,
and in Section 5 we apply the method to data from the thermal pain study. The
paper concludes with a discussion.

##  2 An fMRI Study of Thermal Pain

The data are from an fMRI study of thermal pain; see Woo et al., (2015) for an
in-depth discussion of the data set. In brief, $33$ healthy, right-handed
subjects completed the study (age $27.9\pm 9.0$ years, $22$ females). All
subjects gave informed consent, and the Columbia University Institutional
Review Board approved the study.

Seven runs were administered to each subject during a single session, with
each run consisting of between 58-75 trials. In our analysis the runs were
concatenated and analyzed together. In each trial, thermal stimulations were
delivered to the left inner forearm. Each stimulus lasted $12.5$s, with $3$s
ramp-up and $2$s ramp-down periods and $7.5$s at the target temperature. In
total stimuli were administered at six different temperatures, ranging from
$44.3-49.3^{\circ}$C in increments of $1^{\circ}$C. Each stimulus was followed
by a $4.5-8.5$s pre-rating period, after which subjects rated their intensity
of pain on a scale between $0-100$. Each trial ended with a $5-9$s resting
period.

For each subject, 1845 images were acquired using a 3T Philips Achieva TX
scanner at Columbia University. Structural images were acquired using high-
resolution T1 spoiled gradient recall (SPGR) images. Functional echo planar
images (EPIs) were acquired with repetition time (TR) = $2000$ms, echo time
(TE) = $20$ms, field of view = $224$mm, $64\times 64$ matrix, $3\times 3\times
3$mm3 voxels, $42$ interleaved slices, parallel imaging, and sensitivity
encoding (SENSE) factor $1.5$. For each subject, structural images were co-
registered to the mean functional image using the iterative mutual
information-based algorithm in SPM8111Statistical Parametric Mapping, version
8; http://www.fil.ion.ucl.ac.uk/spm/. Thereafter, images were normalized to
Montreal Neurological Institute (MNI) space using SPM8’s generative segment-
and-normalize algorithm. Prior to preprocessing of functional images, the
first four volumes were removed to allow for image intensity stabilization.
Outliers were identified using the Mahalanobis distance for the matrix of
slice-wise mean and standard deviation values. The functional images were
corrected for differences in slice-timing, and motion corrected using SPM8.
These images were warped to SPM’s normative atlas using warping parameters
estimated from co-registered high-resolution structural images, and smoothed
with an $8$mm full width at half maximum (FWHM) Gaussian kernel. A high-pass
filter of $180$s was applied to the time series data.

For each subject, a voxel-wise general linear model (GLM) analysis was
performed using SPM8. For each temperature, boxcar regressors convolved with
the canonical hemodynamic response function (Lindquist et al.,, 2009), were
constructed to model periods corresponding to the $12.5\,s$ thermal
stimulation and the $11\,s$ rating periods. Other regressors of non-interest
(i.e., nuisance variables) included (a) a run specific intercept; (b) linear
drift across time within each run; (c) the six estimated head movement
parameters (x, y, z, roll, pitch, and yaw), their mean-centered squares, their
derivatives, and squared derivative for each run; (d) indicator vectors for
outlier time points identified based on their multivariate distance from the
other images in the sample; (e) indicator vectors for the first two images in
each run; (f) signals from white matter and ventricle. For each temperature
the estimated regression coefficients corresponding to that temperature
specific regressor at each voxel were combined into a single activation map.
Hence, for each subject we had six different brain maps depicting the
functional response across the brain to each of the temperatures.

In this work, we are primarily focused on the subject-specific mean across the
six difference brain maps to analyze the inter-subject variability. We use
this data to perform functional alignment. In particular, we focus our
analysis on a somatosensory region of the brain; Fig 5 shows  a circular
region of radius 15 voxels across subjects. Within this region one can
identify several local peaks within individual participants that conform to
known somatosensory areas. However, there clearly exists substantial inter-
subject variability in their exact location. The somatosensory region of the
brain was chosen specifically because it is known to be related to pain
(Bushnell et al.,, 1999; Gustin et al.,, 2012; Vierck et al.,, 2013). It is
reliably observed in pain studies as evidenced by meta-analyses of neuroimging
pain studies (see, for example, NeuroSynth (Yarkoni et al.,, 2011)).

##  3 Methods

The goal of image registration is to find a mapping from each voxel in one
image (the floating map) to a corresponding position in another image (the
reference map). To formalize the problem, let
$\mathbf{Y}=(Y(s_{1}),Y(s_{2}),\dots,Y(s_{V}))^{T}$ and
$\mathbf{R}=(R(s_{1}),R(s_{2}),\dots,R(s_{V}))^{T}$ be the floating map and
reference map, respectively, located at the same set of voxel locations
$\mathbf{S}=\\{s_{1},s_{2},\cdots,s_{V}\\}$. This mapping involves operations
(e.g., translation, scaling, rotations, or more complex non-linear
transformations) described by a transformation operator $T$ on the floating
map $\mathbf{Y}$ to align it with the reference $\mathbf{R}$.

Let $\mathbf{w}$ denote the set of parameters characterizing the
transformation function $T:=T(\cdot,\mathbf{w})$. Our goal is to estimate $T$
via $\mathbf{w}$ with proper uncertainty quantification. Hence, we seek a
Bayesian formulation of the registration problem that will allow obtaining
credible intervals for the transformation parameters via posterior MCMC
samples. A conventional approach for Bayesian treatment of such a registration
problem would require likelihood from a generative model that represents the
relationship between the floating and reference maps. One can consider a model
of the form $\mathbf{R}\sim_{d}F(\mathbf{Y}(T),\theta)$, where
$\mathbf{Y}(T)=(Y(T(s_{1})),\ldots,Y(T(s_{V}))^{T}$. This requires the
specification of the full probability distribution $F$, parametrized via
$\theta$. A distribution-free approach is less assumption-driven and more
desirable. Also, such a probability model will typically posit the reference
map $\mathbf{R}$ to be noisier than the floating map. For example, one can
specify a generative model $F$ as

$\mathbf{R}(s)=b\mathbf{Y}(T(s,\mathbf{w}))+\mathbf{e}(s),\;\;\text{for}\;s\in\mathbf{S},$
(1)

where $\mathbf{Y}(T(s,\mathbf{w}))$ is the activation of the floating map at
the transformed location $T(s,\mathbf{w})$, the factor $b$ adjusts for the
discrepancy of intensities between the floating and reference maps, and
$\mathbf{e}(s)$ is the additive noise. It is evident from (1) that
$\mathbf{R}$ is modeled as noisier than $\mathbf{Y}$ due to the errors
$\mathbf{e}$. This is undesirable as in practice the observed floating map is
expected to be equally or more noisy than the reference, and hence generative
models of this form are likely to be misspecified.

Furthermore, the generative model approach becomes problematic when working
with multiple floating maps $\mathbf{Y}_{i}$, as it creates multiple
generative models for the common reference map $\mathbf{R}$. Another possible
model-based solution to the multiple registration problem would be to switch
the roles of $\mathbf{Y}$ and $\mathbf{R}$ in (1), expressing
$\mathbf{Y}(s)=b\mathbf{R}(T_{rv}(s,\mathbf{w}^{\prime}))+$error, where
$T_{rv}:=T_{rv}(\cdot,\mathbf{w}^{\prime})$ denotes the reverse
transformation, mapping $\mathbf{R}$ to $\mathbf{Y}$, parametrized by
$\mathbf{w}^{\prime}$. This formulation easily extends to multiple images
$\mathbf{Y}_{i}$, naturally modeling all the floating maps using the common
reference map $\mathbf{R}$, and mitigates the issue of the reference map being
modeled as more noisy. However, the model fit
$\hat{b}\mathbf{R}(T(s,\hat{\mathbf{w}}))$ removes all subject-specific
features that are not present in the reference. Also this switched approach
does not offer any direct way to register the floating image to the reference.
To remedy this, one often leverages the notion of inverse consistency, jointly
estimating both the forward and the reverse maps $T$ and $T_{rv}$ with the
constraint that their composition is the identity operator. Such an approach
can be seamlessly implemented when using two loss-functions, one
characterizing similarity between $\mathbf{R}$ and $\mathbf{Y}(T)$ and the
other between $\mathbf{Y}$ and $\mathbf{R}(T_{rv})$. However, it is
challenging to make this approach amenable to a fully Bayesian treatment, as
it would involve two different generative models
$\mathbf{Y}|\mathbf{R}(T_{rv})$ and $\mathbf{R}|\mathbf{Y}(T)$ for the same
pair of data $(\mathbf{R},\mathbf{Y})$. As these properties are undesirable,
alternative approaches must be explored to conduct Bayesian registration.

###  3.1 Generalized Bayes framework for registration

For the above mentioned reasons, we abandon a fully model-based formulation
and instead adopt a more flexible general loss function based approach to
measure and correct for relative misalignment of the individual map
$\mathbf{Y}$ with respect to the reference map $\mathbf{R}$. We first consider
a loss-function $\ell(\mathbf{R},\mathbf{Y}(T),\mathbf{w})$ corresponding to
the forward transformation $T$ where
$\mathbf{Y}(T)=(Y(T(s_{1}),\ldots,Y(T(s_{V}))^{T}$ and $\mathbf{w}$ denotes
all the parameters involved in the transformation operator $T$ and the loss
function $\ell$.

If $\pi(\mathbf{w})$ denotes a prior on the parameters, then a posterior
distribution given the loss function $\ell$ will be given by

$\pi(\mathbf{w}\;|\;data)=\frac{\exp\left(-\ell(\mathbf{R},\mathbf{Y}(T),\mathbf{w})\right)\pi(\mathbf{w})}{\int\exp\left(-\ell(\mathbf{R},\mathbf{Y}(T),\mathbf{w})\right)\pi(\mathbf{w})d\mathbf{w}}$
(2)

Posteriors of the form (2) that do not rely on a probability model for the
data but simply a loss function are called Gibbs posteriors or generalized
posteriors. They feature in early works of Vovk, (1990); Shawe-Taylor and
Williamson, (1997); McAllester, (1999); Walker and Hjort, (2001) and their
posterior summaries are referred to as Laplace Type Estimators (LTE) in
Chernozhukov and Hong, (2003) who developed extensive theoretical results on
posterior concentration, asymptotic normality and valid confidence intervals
for such posteriors and estimators. Bissiri et al., (2016) argued that updates
of the form (2) are the only coherent update of a prior belief $\pi$ given the
data and the choice of the loss function $\ell$. Using Gibbs posteriors are
becoming an increasingly common form of Bayesian inference, termed as
generalized Bayes (Grünwald and Mehta,, 2020; Fiksel et al.,, 2021; Rigon et
al.,, 2020).

The loss function approach is compatible with multiple registration to a
common template. When working with multiple floating maps, the total loss is
simply the sum of the losses for the individual maps. For example, registering
$K$ images $\mathbf{Y}_{i},i=1,\ldots,K$ to a common $\mathbf{R}$ entails
simple extension of the loss function to
$\sum_{i=1}^{K}\ell(\mathbf{R},\mathbf{Y}_{i}(T),\mathbf{w}_{i})$ where the
parameters $\mathbf{w}_{i}$ can be image-specific or shared across images,
i.e., $\mathbf{w}_{i}=\mathbf{w}$ depending on the context. Thus the
generalized Bayes approach for registration, by dissociating loss-functions
from likelihoods of generative models, avoids the incompatibility (discussed
in Section 3) of the model-based approach to register multiple floating maps.

We define the loss function using the sum-of-squared differences (SSD) between
the warped image $\mathbf{Y}(T)$ and $\mathbf{R}$,

$\ell(\mathbf{R},\mathbf{Y},b,\phi,\mathbf{w})=\frac{1}{\phi^{2}}\|\mathbf{R}-b\mathbf{Y}(T(\cdot,\mathbf{w}))\|^{2}_{L_{2}}.$
(3)

The generalized-Bayes approach also seamlessly extends to satisfy the property
of inverse-consistency introduced by Christensen and Johnson, (2001) and
Johnson and Christensen, (2002). For this purpose, we tweak the loss function
as follows:

$\begin{array}[]{cc}\ell_{ic}(\mathbf{Y},\mathbf{R},b,b^{\prime},\phi,\mathbf{w},\mathbf{w}^{\prime})=&\ell(\mathbf{R},\mathbf{Y},b,\phi,\mathbf{w})+\ell(\mathbf{Y},\mathbf{R},b^{\prime},\phi,\mathbf{w}^{\prime})+\\\
&\lambda_{rv}\|T_{rv}(T(\cdot))-Id(\cdot)\|_{F}.\end{array}$ (4)

Here $\mathbf{w}^{\prime}$ characterizes the reverse transformation map,
$T_{rv}(\cdot,\mathbf{w}^{\prime})$, from $\mathbf{R}$ to $\mathbf{Y}$ and
$Id(\cdot)$ is the identity map. The scaling parameter $\phi$ controls the
learning rate from the loss function (relative to the prior) in generalized
Bayes and is assumed to be the same for the forward and reverse losses,
thereby assigning equal weight to both. The last quantity in (4) is the
inverse-consistency penalty that constrains the composition of forward and
reverse transformation to be close to the identity map. The hyper-parameter
$\lambda_{rv}$ balances the intensity loss and inverse-consistency penalty.
While there are many potential image similarity measures (Crum et al.,, 2004),
we use SSD as it has been demonstrated as being appropriate for single-modal
brain registration (Simpson et al.,, 2012). However, the general framework
outlined here is agnostic to the choice of loss functions, and other losses
like mutual information loss (Viola and Wells,, 1997; Thevenaz and Unser,,
1998), correlation-based losses (Gruen and Baltsavias,, 1987), can be used
depending on the application.

To implement this approach, decisions need to be made regarding how spatial
interpolation is performed on the discrete images to obtain $\mathbf{Y}(T)$
used by the loss function, choices of transformation operator, reference map,
prior and hyper-parmeters. We discuss these decisions in detail below.

###  3.2 Kriging Interpolation

It is important to note that updating the parameters via optimization, or
sampling using the loss function (3), would require evaluating the warped
images $\mathbf{Y}(T)$ for a given value of the transformation parameters
$\mathbf{w}$. As the activation maps $\mathbf{Y}$ are only observed on a
predefined set of $V$ voxels $\mathcal{V}$, the warped images, $Y(T)$ at the
set of locations $T(\mathcal{V})$, are acquired through spatial interpolation.
While any interpolator can be used for this purpose, we use the Kriging
interpolator (Matheron,, 1963; Banerjee et al.,, 2014; Cressie and Wikle,,
2015)

$\widehat{\mathbf{Y}(T)}=\frac{\mathbf{1}^{T}\mathbf{C}^{-1}\mathbf{Y}}{\mathbf{1}^{T}\mathbf{C}^{-1}\mathbf{1}}\mathbf{1}+\mathbf{K}^{T}\mathbf{C}^{-1}\left(\mathbf{Y}-\frac{\mathbf{1}^{T}\mathbf{C}^{-1}\mathbf{Y}}{\mathbf{1}^{T}\mathbf{C}^{-1}\mathbf{1}}\mathbf{1}\right)$
(5)

where $\mathbf{1}$ is a vector of ones of length $V$ and for a given choice of
a covariance kernel $C$, $\mathbf{C}$ and $\mathbf{K}$ are $V\times V$
matrices with $(i,j)^{t}h$ entries $C(s_{i},s_{j})$ and $C(s_{i},T(s_{j}))$
respectively. We use the popular exponential covariance kernel, i.e.,
$C(s_{i},s_{j})=\sigma^{2}\exp\\{-\rho\|s_{i}-s_{j}\|$ for any $s_{i},s_{j}$,
although any other kernel can be used. The spatial parameters $\sigma$ and
$\rho$ can be estimated within the Bayesian framework jointly with the other
parameters, although increasingly, a priori estimation of these parameters, as
we do here, is being adopted as the pragmatic strategy (Finley et al.,, 2019).

###  3.3 Transformation Operator

The general framework for Bayesian registration proposed via (2) and (4) can
use any family of transformations. We here use the similarity transformation
(Cederberg,, 2001), which consists of translation, scaling, and rotation. For
the case of 2D images this can be expressed as follows:

$T(\mathbf{s})=\mathbf{A}\mathbf{s}+\mathbf{\theta},\mbox{ where
}\mathbf{A}=\begin{pmatrix}\cos(\omega)&-\sin(\omega)\\\
\sin(\omega)&\cos(\omega)\end{pmatrix}\times\begin{pmatrix}\sigma_{x}&0\\\
0&\sigma_{y}\end{pmatrix},$ (6)

where $\mathbf{s}=(s_{x},s_{y})^{T}$,
$\mathbf{\theta}=(\theta_{x},\theta_{y})^{T}\in\mathbb{R}^{2}$,
$\omega\in(-\frac{\pi}{2},\frac{\pi}{2})$, $\sigma_{x}>0$, $\sigma_{y}>0$. The
similarity transformation in (6) is a special case of the general 2D affine
transformation. To see this, we decompose the transformation matrix
$\mathbf{A}$ as

$\displaystyle\mathbf{A}=\begin{pmatrix}A_{11}&A_{12}\\\
A_{21}&A_{22}\end{pmatrix}$ (7)
$\displaystyle=\begin{pmatrix}\cos(\omega)&-\sin(\omega)\\\
\sin(\omega)&\cos(\omega)\end{pmatrix}\begin{pmatrix}1&\frac{A_{12}\cos(\omega)+A_{22}\sin(\omega)}{A_{22}\cos(\omega)-A_{12}\sin(\omega)}\\\
0&1\end{pmatrix}\begin{pmatrix}\text{sign}(A_{11})\sqrt{A_{11}^{2}+A_{21}^{2}}&0\\\
0&A_{22}\cos(\omega)-A_{12}\sin(\omega)\end{pmatrix},$

where $\omega=\arctan(A_{21}/A_{11})$. We note that the expression above is
the result of the QR decomposition. The first term is the rotation matrix. It
is followed by the shearing matrix. The last term is the scaling matrix. The
equivalence between Eqns. (6) and (7) can now be derived by adding constraints
to eliminate the shearing effect, i.e., setting
$\frac{A_{12}\cos(\omega)+A_{22}\sin(\omega)}{A_{22}\cos(\omega)-A_{12}\sin(\omega)}=0$,
or equivalently $A_{11}A_{12}+A_{21}A_{22}=0$. Under this constraint, the
decomposition can be simplified as

$\mathbf{A}=\begin{pmatrix}A_{11}&A_{12}\\\
A_{21}&A_{22}\end{pmatrix}=\begin{pmatrix}\cos(\omega)&-\sin(\omega)\\\
\sin(\omega)&\cos(\omega)\end{pmatrix}\begin{pmatrix}\text{sign}(A_{11})\sqrt{A_{11}^{2}+A_{21}^{2}}&0\\\
0&A_{22}/\cos(\omega)\end{pmatrix}.$ (8)

Therefore, there exists a bijection mapping between
$(\sigma_{x},\sigma_{y},\omega)^{T}$ and $\text{vec}(A)$ under the constraint
that $A_{11}A_{12}+A_{21}A_{22}=0$. Furthermore, if we are interested in the
positive scalings, $\sigma_{x}>0$, $\sigma_{y}>0$, this can be achieved by
setting $A_{11}>0$ and $A_{22}>0$. We will exploit this connection to
formulate the prior for the transformation parameters in the next Section.

###  3.4 Regularization Prior

There are several obstacles in optimizing for transformation parameters in the
general registration problem. First, it is an inherently ill-posed problem. A
classic example where two distinct registration routines might lead to the
same registration result was discussed in Fischer and Modersitzki, (2008),
where they showed a subtle example that either translation or rotation can
yield the same optima. Second, it is difficult to avoid unwanted local optima.
Both of these difficulties can be addressed by including a regularization, or
penalty, term. In this Section we adapt this idea by introducing
regularization priors for our generalized Bayes framework.

The generalized (Gibbs) posterior distribution of the warping parameters takes
the form:

$\displaystyle
p(T,T_{rv},\phi,b,b^{\prime}|\mathbf{Y},\mathbf{R})\propto\exp{\left\\{-\ell_{ic}(\mathbf{Y},\mathbf{R},b,b^{\prime},\phi,\mathbf{w},\mathbf{w}^{\prime})\right\\}}\pi(T,T_{rv},\phi,b,b^{\prime}),$
(9)

where $\pi(T,T_{rv},\phi,b,b^{\prime})$ is the regularization prior. We
specify $\pi(T,T_{rv},\phi,b,b^{\prime})=$
$\pi(T|\phi)\pi(T_{rv}|\phi)\pi(b^{\prime}|\phi)\pi(b|\phi)\pi(\phi)$ for
choices of $\pi(T|\phi)$, $\pi(T_{rv}|\phi)$, $\pi(b|\phi)$ and
$\pi(b^{\prime}|\phi)$ which leads to the following form of the negative log-
posterior,

$\displaystyle\frac{1}{2\phi^{2}}\|\mathbf{R}-b\mathbf{Y}(T)\|_{L_{2}}^{2}+\frac{1}{2\phi^{2}}\|\mathbf{Y}-b\mathbf{R}(T_{rv})\|_{L_{2}}^{2}+\lambda_{rv}\|T_{rv}(T(\cdot))-Id(\cdot)\|_{F}$
(10)
$\displaystyle+\frac{\lambda_{T}}{\phi^{2}}\text{Reg}(T)+\frac{\lambda_{b}}{\phi^{2}}(\log(b)-\log(b_{0}))^{2}+\frac{\lambda_{T_{rv}}}{\phi^{2}}\text{Reg}(T_{rv})+\frac{\lambda_{b^{\prime}}}{\phi^{2}}(\log(b^{\prime})-\log(b^{\prime}_{0}))^{2},$

where the regularizer Reg can be chosen in different ways. The terms
$\lambda_{b}$ and $\lambda_{T}$ (and analogously, $\lambda_{b^{\prime}}$ and
$\lambda_{T_{rv}}$) are the tuning parameters, which control the trade-off
between the fitting of the registration likelihood and the regularization
priors. The factor $(\log(b)-\log(b_{0}))^{2}$ in the regularization part
regularizes the intensity-correction term, and the term $b_{0}$ is one of the
hyper-parameters whose choice will be described in the next Section. To gain
intuition on the necessity of this term, without penalizing values of $b$ away
from some $b_{0}$, the MCMC may become stuck in loss function valleys centered
around a local minima. For example, around $b=0$, the loss function will be
flat for all values of the transformation parameters $w$. The prior for $b$
leading to (10) is given by $\log(b)|\phi\sim
N(\log(b_{0}),~{}\frac{\phi^{2}}{\lambda_{b}})$. The regularization of $T$,
$\text{Reg}(T)$, is taken to measure the distance between the coordinates
transformed by the proposed $T$ and the some initial warping map $T_{0}$
defined by hyper-parameters $A_{0}$ and $\mathbf{\theta}_{0}$. The selection
of these hyper-parameters is critical to ensure identifiability of the
registration, and is discussed in detail in Section 3.5. Formally, the
regularizer of $T$ is represented in the following way:

$\text{Reg}(T)=\sum_{\mathbf{s}\in\mathcal{R}}\|T(\mathbf{s})-T_{0}(\mathbf{s})\|_{L_{2}}^{2}=\sum_{\mathbf{s}\in\mathcal{R}}\|(A-A_{0})\mathbf{s}+\mathbf{\theta}-\mathbf{\theta}_{0}\|^{2}_{L_{2}}.$
(11)

The expression of the regularizer can be further simplified as:

$\displaystyle\text{Reg}(T)$
$\displaystyle=\text{trace}\\{[\mathbf{s}_{x}~{}\mathbf{s}_{y}~{}\mathbf{1}](M-M_{0})(M-M_{0})^{T}[\mathbf{s}_{x}~{}\mathbf{s}_{y}~{}\mathbf{1}]^{T}\\}$
$\displaystyle=\text{trace}\\{(M-M_{0})^{T}\Sigma_{s}(M-M_{0})\\},$

where $M=[A~{}\mathbf{\theta}]^{T}=\begin{pmatrix}A_{11}&A_{12}&\theta_{x}\\\
A_{21}&A_{22}&\theta_{y}\end{pmatrix}^{T}$,
$M_{0}=[A_{0}~{}\mathbf{\theta}_{0}]^{T}$, and $\Sigma_{s}=$
$\begin{pmatrix}\mathbf{s}_{x}^{T}\mathbf{s}_{x}&\mathbf{s}_{x}^{T}\mathbf{s}_{y}&\mathbf{s}_{x}^{T}\mathbf{1}\\\
\mathbf{s}_{x}^{T}\mathbf{s}_{y}&\mathbf{s}_{y}^{T}\mathbf{s}_{y}&\mathbf{s}_{y}^{T}\mathbf{1}\\\
\mathbf{s}_{x}^{T}\mathbf{1}&\mathbf{s}_{y}^{T}\mathbf{1}&\mathbf{1}^{T}\mathbf{1}\\\
\end{pmatrix}$.

Hence, using the prior $\text{vec}(M)|\phi\sim
N(\text{vec}(M_{0}),~{}\frac{\phi^{2}}{\lambda_{T}}\mathbf{I}\otimes\Sigma_{s}^{-1})$,
where $\otimes$ denotes the Kronecker product and
$\mathbf{I}\in\mathbb{R}^{2\times 2}$ the identity matrix, ensures that the
Gibbs posterior is of the form (10). Note that, $M$ has $6$ parameters and our
transformation is only defined by $5$ parameters
$\omega,\sigma_{x},\sigma_{y},\theta_{x},\theta_{y}$. As discussed in Section
3.3, there is a one-to-one relationship between our parametrizations and the
general 2D affine transformation parametrized by $M$ under the constraint
$A_{11}A_{12}+A_{21}A_{22}=0$. Hence, we only use the multivariate normal
prior for the 3-dimensional subset of $A$ along with $\mathbf{\theta}$ to
correspond to the priors on the rotation, scaling and translation parameters.

Similarly, we define the regularizer of $T_{rv}$ as

$\displaystyle\text{Reg}(T_{rv})=\text{trace}\\{(M^{\prime}-M^{\prime}_{0})^{T}\Sigma_{s}(M^{\prime}-M^{\prime}_{0})\\},$
(12)

where $M^{\prime}=[A^{\prime}~{}\mathbf{\theta}^{\prime}]^{T}$,
$T_{rv}(s)=A^{\prime}\mathbf{s}+\theta^{\prime}$.
$M^{\prime}_{0}=[A^{\prime}_{0}~{}\mathbf{\theta}^{\prime}_{0}]^{T}$ is the
prior transformation parameters of $T_{rv}$. To meet the requirement of
inverse-consistency, the prior mapping, $T_{rv0}$ is constrained to the
inverse of $T_{0}$. In other words, we have
$A^{\prime}_{0}(A_{0}\mathbf{s}+\theta_{0})+\theta^{\prime}_{0}=\mathbf{s}$.

Finally, we use the standard reference prior for the scale parameter $\phi$,
i.e., $\pi(\phi)\propto 1/\phi^{2}$. Combining all the pieces in the
generalized Bayes framework, the Gibbs posterior distribution of the warping
parameters has the following form,

$\displaystyle p(\mathbf{\theta},\mathbf{\sigma},\omega)\propto$ (13)
$\displaystyle\phi^{-2V-14}\exp\Bigg{\\{}-\frac{1}{2\phi^{2}}\|\mathbf{R}-b\mathbf{Y}(T)\|_{L_{2}}^{2}-\frac{1}{2\phi^{2}}\|\mathbf{Y}-b^{\prime}\mathbf{R}(T_{ic})\|_{L_{2}}^{2}$
$\displaystyle-\frac{\lambda_{b}}{2\phi^{2}}((\log(b)-\log(b_{0}))^{2})-\frac{\lambda_{T}}{2\phi^{2}}\text{trace}\\{(M-M_{0})^{T}\Sigma_{s}(M-M_{0})\\}$
$\displaystyle-\frac{\lambda_{b^{\prime}}}{2\phi^{2}}((\log(b^{\prime})-\log(b^{\prime}_{0}))^{2})-\frac{\lambda_{T_{rv}}}{2\phi^{2}}\text{trace}\\{(M^{\prime}-M^{\prime}_{0})^{T}\Sigma_{s}(M^{\prime}-M^{\prime}_{0})\\}\Bigg{\\}}.$

Notice that we have dropped the explicit penalty for inverse-consistency, as
it is embedded in the choice of prior transformations.

###  3.5 Corresponding Feature-based Prior Estimation

Informed selection of the regularization hyper-parameters $b_{0}$ and $M_{0}$
is critical to ensure identifiability of the registration problem. In this
Section we outline an algorithm for prior (hyper-parameter) estimation based
on corresponding features. The loss function (3) is solely based on the
intensity information of the images, and does not use the landmark information
of activation maps which is an alternate branch in the image registration. The
prior (hyper-parameters) of the transformation can be estimated by
incorporating information about the landmarks. This approach does not
encounter computational burden with a moderate number of landmarks. In this
work, we define the landmarks as the regional peaks (maxima) in clusters
consisting of eight-adjacent voxels.

Let $\mathbf{p}^{R}=\\{s_{i}^{R}\in\mathbb{R}^{2}|i=1,\cdots,N_{R}\\}$ and
$\mathbf{p}^{Y}=\\{s_{i}^{Y}\in\mathbb{R}^{2}|i=1,\cdots,N_{Y}\\}$ be the
landmarks of reference map and the floating map, respectively. $N_{R}$ and
$N_{Y}$ are the number of landmarks in the reference map and floating map,
respectively. As a first step, the coarse sub-region of interest of the
reference map is determined manually based on scientific beliefs. This sub-
region contains a subset of landmarks,
$\hat{\mathbf{p}}^{R}\in\mathbf{p}^{R}$, of size denoted by $\hat{N}_{R}$. We
aim to find the corresponding subset of the landmarks $\hat{\mathbf{p}}^{Y}$
from $\mathbf{p}^{Y}$. Many conventional methods, including Geometric Hashing
(Mian et al.,, 2006) and Iterative Closest Point (Yang and Medioni,, 1992;
Besl and McKay,, 1992), work for point set matching. However, these methods
fail to capture the intensity features of the maps. In our application, as the
size of the set of landmarks are relatively small (i.e., $N_{Y}<20$ and
$\hat{N}_{R}<5$), it is computationally efficient to apply the Brute Force
Search method that tries all possible $\hat{N}_{R}$ dimensional subsets of
$\mathbf{p}^{Y}$ to determine the matching landmarks.

The similarity transformation can be estimated using Procrustes analysis. Here
we apply a modified version of Procrustes analysis, namely ABC Procrustean
Algorithm (Awange et al.,, 2008). This version is preferable to the ordinary
Procrustes algorithm as it is capable of finding the distinct scale parameters
in each axial direction as opposed to the estimation of uniform scales.

The objective function used to search for correspondence is selected as the
sum of squared distance between the query landmarks and the transformed
candidates. Again, we need to add a constraint term to address the ill-posed
registration problem, as discussed in the beginning of this Section. Formally,
the mechanism can be formulated as follows,

$\displaystyle||\hat{\mathbf{p}}^{R}-T_{\alpha}(\hat{\mathbf{p}}^{Y})||^{2}+||\hat{\mathbf{p}}^{Y}-T^{-1}_{\alpha}(\hat{\mathbf{p}}^{R})||^{2}\leq
2d$ (14) subject to $\displaystyle
S(\hat{\mathbf{p}}^{Y}-T_{\alpha}(\hat{\mathbf{p}}^{Y}))+S(\hat{\mathbf{p}}^{R}-T^{-1}_{\alpha}(\hat{\mathbf{p}}^{R}))<\alpha,$

using the transformation operator
$T_{\alpha}(\mathbf{s})=\mathbf{A}_{\alpha}\mathbf{s}+\mathbf{\theta}_{\alpha}$.

The fact that the activation maps are discrete realizations of continuous maps
motivates the use of the distance less than the threshold $d$, instead of
minimizing the distance between landmarks. The discreteness of a map results
in the bias of the coordinates of landmarks, as the underlying location of
landmarks can be anywhere within a pixel. We suggest to take $d=\hat{N}_{R}$
which relaxes the landmark within a voxel. The objective function creates a
set of candidates of correspondence, and we choose the optimal correspondence
based on the intensity loss.

The constraint $S$ is motivated by diffusion registration (Fischer and
Modersitzki,, 2003) and has been applied to affine models by Chumchob and
Chen, (2009). By defining a function
$D\colon\mathbb{R}^{2}\rightarrow\mathbb{R}^{2}$ with
$D(\mathbf{p})=(D_{1}(\mathbf{p}),D_{2}(\mathbf{p}))^{T}=\mathbf{p}-T(\mathbf{p})$,
it can be expressed as follows:

$\displaystyle S(\mathbf{p}-T(\mathbf{p}))=\frac{1}{2}\int(|\nabla
D_{1}(\mathbf{p})|^{2}+|\nabla
D_{2}(\mathbf{p})|^{2})d\mathbf{p}=\|\mathbf{I}-\mathbf{A}||_{F}^{2},$ (15)

where $\|\cdot\|_{F}$ denotes the Frobenius norm, $\mathbf{I}$ is the identity
matrix, $\mathbf{A}$ is defined in Eq. (6). The other application of the
constraint defined in Eq. (15) is to determine the orientation of
$\hat{\mathbf{p}}^{Y}$. For example, in the case of $\hat{N}^{R}=3$,
$\hat{\mathbf{p}}^{Y}=\\{s_{1}^{Y},s_{2}^{Y},s_{3}^{Y}\\}$, there are six
possible correspondences to $\hat{\mathbf{p}}^{R}$, and the transformation can
be derived for each of them. The mismatched correspondence may result in
excessive transformation, which can cause redundant rotation and scaling.
Therefore, the one minimizing
$S(\hat{\mathbf{p}}^{Y}-T_{\alpha}(\hat{\mathbf{p}}^{Y}))$ is used to
determine the orientation of $\hat{\mathbf{p}}^{Y}$.

The threshold $\alpha$ in Eq. 14 is chosen automatically using the Photometric
Error Criterion (PEC) (Brunet et al.,, 2010). For each $\alpha$, the
optimization problem in Eq. (14) is used to determine the corresponding
features as well as the associated transformation. To evaluate the hyper-
parameter $\alpha$, the criterion takes advantage of intensity information,
which is defined as:

$\mathcal{C}(\alpha)=\frac{1}{|\mathcal{G}|}\sum_{\mathbf{p}\in\mathcal{G}}||\mathbf{R}(\mathbf{p})-b\mathbf{Y}(T_{\alpha}(\mathbf{p}))||^{2}+\frac{1}{|\mathcal{G^{\prime}}|}\sum_{\mathbf{p}\in\mathcal{G^{\prime}}}||\mathbf{Y}(\mathbf{p})-b^{\prime}\mathbf{R}(T^{-1}_{\alpha}(\mathbf{p}))||^{2}.$
(16)

Here $\mathcal{G}$ and $\mathcal{G^{\prime}}$ are the regions of interest in
$\mathbf{R}$ and $\mathbf{Y}$ respectively, and $|\cdot|$ defines the size of
the region. The discrepancy of scale between the target map and the reference
map is corrected by a scaling factor $b$. Note that $b$ can be estimated as
the regression coefficient of $\mathbf{R}(\mathbf{p})$ against
$\mathbf{Y}(T_{\alpha}(\mathbf{p}))$ without an intercept as the criterion in
Eq. (16) is denoted as the mean squared error of the regression.

Eq. 16 is the summation of mean squared error of regressing $\mathbf{R}$
against $\mathbf{Y}(T_{\alpha})$ and $\mathbf{Y}$ against
$\mathbf{R}(T^{-1}_{\alpha})$, where $b^{\prime}$ is the regression
coefficient. This estimate $b_{0}$ and $b^{\prime}_{0}$ is used as a hyper-
parameter in the prior of $b$ and $b^{\prime}$ in Section 3.5.

We may further define the robust photometric error criterion (RPEC) by making
a heavy-tailed distribution assumption on the residuals, for instance, the
Student T distribution with a small degree freedom. The selected $\alpha$
minimizes the criterion in Eq. (16), and produces the feature correspondence
and the estimation of transformation $M_{0}$ and $M^{\prime}_{0}$ which will
be considered as the hyper-parameter in the prior for $M$ and $M^{\prime}$
respectively in Section . We summarize the full steps of prior estimation in
Algorithm (1).

###  3.6 Choice of Reference Map

A critical decision when performing registration is deciding on an appropriate
reference map. One option is to use the group mean map across all floating
maps. However, this will generally make the reference map blurry and the
estimated transformation largely variant. In our application, the reference
subject is chosen based on scientific knowledge. For example, researchers may
choose the functional map that best represents the areas of interest to be the
reference.

It is worth noting that instead of aligning the whole brain functional map, we
instead focus on estimating the transformation with respect to local regions
of interest. In this way, estimation will not be unduly influenced by other
nuisance regions that do not contain signal of interest. The query map cropped
by a rectangle bounding box from the reference map describes the local area of
interest; see Figure (1) for an example. The transformed bounding boxes crop
the floating maps to match and register to the query map. As the affine
transformation we focus on is linear, the variability of bounding boxes on
floating maps is equivalent to the variability of registration.

###  3.7 Model Fitting

The inference on the parameters is based on the draws from the posterior
distribution. We use the Markov Chain Monte Carlo (MCMC) algorithm offered in
the Stan software with default settings to simulate from the posterior
distributions. Stan (Carpenter et al.,, 2017) is implemented with the No-U-
Turn sampler (NUTS) (Hoffman and Gelman,, 2014). We refer readers interested
in details of the algorithm to the original work. Three distinct chains are
fitted and initialized at values close to the prior estimates in order to
accelerate the model fitting.

The only hyper-parameters that are not determined in Section 3.5 are
$\lambda_{b}$, $\lambda_{T}$, $\lambda_{b^{\prime}}$ and $\lambda_{T_{rv}}$,
which control the amount of regularization. It is common to choose such hyper-
parameters using cross-validation. However the use of cross-validation
embedded into MCMC runs is computationally intensive, and it requires a large
amount of parallel computation to cover all possible splits of the data.
Vehtari et al., (2017) proposed an approach for approximating leave-one-out
cross-validation using Pareto smoothed importance sampling (PSIS-LOO) which
can be easily computed. Compared to the Watanabe-Akaike information (WAIC;
Watanabe, (2010)), we prefer PSIS-LOO as it is more robust in the finite case
with weak priors or influential observations.

Along with the PSIS, we also calculate the potential scale reduction statistic
(Gelman and Rubin,, 1992), $\hat{R}$, as a metric of sampling convergence.
Starting with a relatively high value of the hyper-parameters, we gradually
fit the model with decreasing $\lambda$ until $\hat{R}<1.01$. This step keeps
the model away from unwanted local modes. For each choice of ($\lambda_{b}$,
$\lambda_{T}$ $\lambda_{b^{\prime}}$, $\lambda_{T_{rv}}$) we compute PSIS and
select the one minimizing PSIS.

The main source of computational complexity is from the Kriging interpolation,
described in Section (3.2). The computational complexity is in the order of
$O(V^{3})$. While it is computationally manageable in our experiment due to
the relative small $V$, the computation may be intractable for large $V$. In
the case of large $V$, we recommend to apply the nearest-neighbor
approximation introduced by Datta et al., (2016), which will efficiently
approximate the Kriging interpolation in the order of $O(Vm^{3})$ for
neighborhood size of $m<<V$. As we are independently registering each
subject’s map, it is completely parallelizable across subjects.

The models were fit in parallel across the choice of regularization
parameters. Computations were carried out on a high-performance computing
(HPC) cluster with an AMD EPYC 7702 64-core processor at 2.0GHz. We ran 10000
sampling iterations and set the initial 2000 iterations as the burn-in stage.
The average computing time was 10 minutes for the burn-in stage and 20 minutes
for the sampling stage. STAN has in-built diagnostics (e.g., the Gelman-Rubin
potential scale reduction factor) that were used for determining convergence.
A software implementation of our method can be found at
https://github.com/gqwang1001/BayesianFunctionalRegistrationOfFMRIMaps.

##  4 Simulation

To illustrate the validity of our method we perform a simulation study using
maps synthesized from the real data. The left panel of Figure 1A shows an
example functional activation map. The query part of the map is manually
cropped and shown within the dashed rectangle and in panel B. The floating
map, shown in panel C, is generated by warping the reference map with the
similarity transformation (including translation, scaling and rotation). The
underlying parameters are translations $\theta=(2,-5)$, scalings
$\sigma=(0.8,1.2)$, and rotation $\omega=\pi/12$. White noise,
$10^{-5}\epsilon$ with $\epsilon\sim N(0,1)$, is subsequently added to the
map.

We compare the learned prior information with the posterior in Figure 2A,
where the region surrounded by the dashed box represents the prior information
and the solid box crops the matched region transformed by the posterior
parameters. The boxes are generated by transforming the border box of the
reference map using the prior estimates and median of posterior samples,
respectively, of the transformation parameters.

The models were fit using a sequence of regularization parameters
($\lambda_{b}$, $\lambda_{T}$, $\lambda_{b^{\prime}}$, $\lambda_{T_{rv}}$) in
Eq. (9). The expected log predictive density (ELPD) was estimated on each
model. Figure 2B plots the -2 $\times$ ELPD which makes it on the deviance
scale. The samples produced by the model with the lowest -2$\times$ ELPD were
used to estimate the posterior distribution. The uncertainty of transformation
can be demonstrated as posterior distribution of the translation, scaling and
rotation parameters which are shown on Figure 3 with the truth superimposed as
the dashed vertical line, where we can see the posterior samples successfully
recover the truth. We further ran 100 simulations with different underlying
transformation. Figure 4 shows the plot of posterior means against the true
transformation. We could see the method recovered the truth successfully in
each case.

To explore the models performance under a mis-specified model, the map is
further deformed nonlinearly by elastic distortions, following the steps
described in Simard et al., (2003). First, the deformation field, $\Delta x$
and $\Delta y$ on $x$ and $y$ directions respectively, are generated from the
uniform distribution, $\text{Unif}(-1,1)$, at each voxel location. We then
convolve the deformation field with the Gaussian filter with standard
deviation $\sigma_{d}=4$. It is subsequently normalized and multiplied by a
scaling factor $\alpha_{d}=60$. Figure S1 shows the posterior mean plotted
against the true affine parameters. Similar to the results shown in Figure 4,
the method successfully recovers the truth, while the deviance is larger due
to the presence of nonlinear generalization.

##  5 Real Data Results

We use the data described in Section 2 to perform functional alignment. Figure
5 shows functional activation maps of the region of interest across the $33$
subjects. Note that within this region there are multiple peaks. However,
there also exists substantial inter-subject variability in their exact
location. For example, comparing subjects 14 and 20 suggests a potential
transformation consisting of rotation and shifting. Each of these functional
activation maps are processed through the adaptive thresholding introduced by
Bradley and Roth, (2007). The landmarks are then identified as the regional
maxima of the neighborhood of 8 adjacent locations.

The same reference map is utilized as described in the simulation; see Figure
1. This particular map was chosen as it was determined to best represent the
area of interest. The subject specific maps are shown in Figure 6, and the
estimated corresponding areas are indicated using the dashed boxes. Figure 7
shows the matching area warped by the posterior mean of transformation, along
with the 95% credible region of the transformations are superimposed in the
gray-shaded region. The credible region of transformation is found using the
density-based clustering algorithm (DBSCAN), introduced by Ester et al.,
(1996). The DBSCAN algorithm is the unsupervised clustering algorithm that
separates the clusters based on the density. We refer readers to the original
works for further details. To apply DBSCAN, we specify the epsilon region with
at least 5 samples and a grid search of the size of the epsilon neighborhood
to assure the only cluster captures about 95% of the posterior samples. The
shaded region in the figure is then obtained by warping the searching box with
each posterior samples of the transformation. The inter-subject variability in
registration can be visualized using the estimated credible regions. For
instance, subject 21 exhibits a high degree of uncertainty, while subject 14
exhibits relatively low uncertainty. For subject 6, the estimated area falls
outside the activation map. This is caused by the fact that we are analyzing a
subset of the data corresponding to a circular ‘searchlight’. For this subject
the best corresponding features lies on the edge of the searchlight, leading
to a situation where the corresponding subregion falls outside of the
searchlight range.

We contrast the results to those obtained using a standard non-linear
deformation-based method. We applied the pipeline implemented in the image
processing toolbox in MATLAB (version R2020b). Affine registration was first
performed between each subject’s activation map and the reference map,
followed by diffeomorphic demons registration (Vercauteren et al.,, 2009).
Figure 8 demonstrates the warping results. Seven out of thirty-three subjects,
including subject 3, 5, 6, 12, 19, 25, 27, failed to register to the reference
map. This is due to the diffuse activation maps. In addition, the standard
approach was unable to provide the types of uncertainty measures of the
registration that our approach provides.

To evaluate the performance of the functional alignment procedure, we
performed a second-level (across-subjects) analysis to test for significant
temperature related effects. Figure 9 shows results comparing our approach to
the standard approach, which first spatially registers images to a common
template and then identifies activated regions, and the diffeomorphic demons
approach. Note that all approaches perform an initial anatomical registration
into a common anatomical space and then compute activation maps. The
difference lies in that our proposed approach and the demons approach
thereafter perform a secondary functional registration of the derived
activation maps. Studying the voxel-wise group means and standard deviations,
as well as t-statistics show increased similarity in the activation profiles
across subjects. After functional alignment, three peaks corresponding to
posterior insula, SII, and opercular are readily apparent. In addition, the
variation is decreased and more focused on the three regions-of-interest. This
indicates that while the locations of the peaks are now consistent across
subjects, the inter-subject differences in intensity remain. This is important
for testing for brain-behavior correlations across subjects. Finally, the
t-statistic is higher, indicating increased sensitivity for group-level
inference. The results obtained using the demons algorithm show a similar
pattern as seen using our approach. However, the t-statistic is significantly
lower, indicating a lower sensitivity for group-level inference. Together, the
results show the benefits of our functional registration approach.

##  6 Discussion

In this work, we propose a framework for reducing misalignment across
individuals in functional brain systems. Our approach has the ability to
register local brain features, defined as local peaks in the activation maps,
providing a means to align functional features prior to inter-subject
statistical analysis. Our generalized Bayes framework combines intensity-based
information through the loss function and the feature-based information
through the prior estimation, allowing us to conduct inference on the
transformation via the posterior samples. In an application to fMRI data from
a study of thermal pain, we found that the registration contributes to
increased sensitivity for group-level inference.

In its current implementation our approach has some limitations. First, in the
case that significant non-linear differences between activation maps change
the relative position between landmarks, our proposed method will likely fail.
Using a more general class of transformations could address this issue.
Second, the choice of cropping boxes can affect the performance of
registration. For example, in cases where the cropping box covers redundant
area (e.g., including extra peak activation voxels that does not match the
floating maps), the performance of the proposed method may suffer due to the
outlier region. Third, we currently choose the reference image subjectively or
based on the results of a meta-analysis. Therefore, in order to apply the
method researchers need to check to ensure that the chosen reference pattern
is shared across subjects. The proposed method can encounter difficult
situations when the reference map is poorly defined. In future work, we will
extend our approach to model the target image as a latent map to be estimated
in a data-driven manner.

In this paper we have presented a 2-dimensional implementation of our proposed
approach. However, we note that a 3-dimensional version could be constructed
in an analogous manner. We began with the 2D version for illustration
purposes, as it allows us to present a simplified description of the approach.
In practice when applying the proposed approach to 3D data we would recommend
a 3D implementation.

Our approach uses functional registration of activation maps to account for
residual differences in brain anatomy after spatial registration to MNI space
has been performed. Once these functional transformations have been derived,
it is possible to apply these transformations to the original 4D fMRI data and
use these to compute functional connectivity. In addition, if the
transformations are capturing true local differences in functional topology
these transformations should carry over to other datasets collected on the
same subject. For example, one could use activation maps from a task-fMRI to
derive local transformations and then apply these same transformation to
resting-state data, which in turn can be used to compute functional
connectivity. Another potential application of our proposed approach is
predictive modeling. Here activation maps are predictors used to model some
behavioral or clinical outcome. Functional registration promises to align
features in the model which can help improve prediction accuracy. In future
work we will seek to explore how functional registration can help improve
population-level models of the functional brain representations that can be
used to predict behavior, performance, clinical status and prognosis. We
hypothesize that due to functional misalignment, features used in predictive
models may not be properly aligned, negatively impacting predictive accuracy.
Using our approach as a preliminary feature alignment step could potentially
rectify this shortcoming.

Furthermore, as our data experiment measures multiple brain maps at six
different temperatures, in future work we will model the within-subject
variability using the registration and explore its association with other
effects such as temperature. This will allow us to explore how the spatial
extent of brain activation changes as a function of stimulus intensity. This
provides a means to move past simply looking at spatial intensity in a
specific voxel as a measure of brain activation. Another possible extension is
evaluating the multilevel variability through the inference on the within-
subject and between-subject transformation.

## Acknowledgments

The work presented in this paper was supported in part by NIH grants R01
EB016061 and R01 EB026549 from the National Institute of Biomedical Imaging
and Bioengineering and R01 MH116026 from the National Institute of Mental
Health.

## References

  * Allison et al., (1999) Allison, T., Puce, A., Spencer, D. D., and McCarthy, G. (1999).  Electrophysiological studies of human face perception. i: Potentials generated in occipitotemporal cortex by face and non-face stimuli.  Cerebral cortex, 9(5):415–430.
  * Amunts et al., (2000) Amunts, K., Malikovic, A., Mohlberg, H., Schormann, T., and Zilles, K. (2000).  Brodmann’s areas 17 and 18 brought into stereotaxic space—where and how variable?  Neuroimage, 11(1):66–84.
  * Awange et al., (2008) Awange, J. L., Bae, K. H., and Claessens, S. J. (2008).  Procrustean solution of the 9-parameter transformation problem.  Earth, Planets and Space, 60(6):529–537.
  * Banerjee et al., (2014) Banerjee, S., Carlin, B. P., and Gelfand, A. E. (2014).  Hierarchical Modeling and Analysis for Spatial Data, Second Edition.  CRC Press, Philadelphia, PA.  ID: 1598075.
  * Besl and McKay, (1992) Besl, P. J. and McKay, N. D. (1992).  A Method for Registration of 3-D Shapes.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 14(2):239–256.
  * Bissiri et al., (2016) Bissiri, P. G., Holmes, C. C., and Walker, S. G. (2016).  A general framework for updating belief distributions.  Journal of the Royal Statistical Society. Series B, Statistical methodology, 78(5):1103.
  * Bradley and Roth, (2007) Bradley, D. and Roth, G. (2007).  Adaptive Thresholding using the Integral Image.  Journal of Graphics Tools, 12(2):13–21.
  * Brunet et al., (2010) Brunet, F., Bartoli, A., Navab, N., and Malgouyres, R. (2010).  Pixel-based hyperparameter selection for feature-based image registration.  VMV 2010 - Vision, Modeling and Visualization, pages 33–40.
  * Bushnell et al., (1999) Bushnell, M., Duncan, G., Hofbauer, R., Ha, B., Chen, J.-I., and Carrier, B. (1999).  Pain perception: is there a role for primary somatosensory cortex?  Proceedings of the National Academy of Sciences, 96(14):7705–7709.
  * Carpenter et al., (2017) Carpenter, B., Lee, D., Brubaker, M. A., Riddell, A., Gelman, A., Goodrich, B., Guo, J., Hoffman, M., Betancourt, M., and Li, P. (2017).  Stan: A probabilistic programming language.  Journal of Statistical Software, 76(1):1–32.
  * Cederberg, (2001) Cederberg, J. N. (2001).  A Course in Modern Geometries.  Undergraduate Texts in Mathematics. Springer New York, New York, NY.
  * Chen et al., (2014) Chen, P. H., Guntupalli, J. S., Haxby, J. V., and Ramadge, P. J. (2014).  Joint SVD-Hyperalignment for multi-subject FMRI data alignment.  IEEE International Workshop on Machine Learning for Signal Processing, MLSP, pages 1–6.
  * Chernozhukov and Hong, (2003) Chernozhukov, V. and Hong, H. (2003).  An mcmc approach to classical estimation.  Journal of Econometrics, 115(2):293–346.
  * Christensen and Johnson, (2001) Christensen, G. and Johnson, H. (2001).  Consistent image registration.  IEEE Transactions on Medical Imaging, 20(7):568–582.
  * Chumchob and Chen, (2009) Chumchob, N. and Chen, K. (2009).  A robust affine image registration method.  International Journal of Numerical Analysis and Modeling, 6(2):311–334.
  * Conroy et al., (2013) Conroy, B. R., Singer, B. D., Guntupalli, J. S., Ramadge, P. J., and Haxby, J. V. (2013).  Inter-subject alignment of human cortical anatomy using functional connectivity.  NeuroImage, 81:400–411.
  * Cressie and Wikle, (2015) Cressie, N. and Wikle, C. K. (2015).  Statistics for spatio-temporal data.  John Wiley & Sons.
  * Crum et al., (2004) Crum, W. R., Hartkens, T., and Hill, D. L. (2004).  Non-rigid image registration: Theory and practice.  British Journal of Radiology, 77(SPEC. ISS. 2).
  * Datta et al., (2016) Datta, A., Banerjee, S., Finley, A. O., and Gelfand, A. E. (2016).  Hierarchical nearest-neighbor gaussian process models for large geostatistical datasets.  Journal of the American Statistical Association, 111(514):800–812.
  * Duncan et al., (2009) Duncan, K. J., Pattamadilok, C., Knierim, I., and Devlin, J. T. (2009).  Consistency and variability in functional localisers.  Neuroimage, 46(4):1018–1026.
  * Ester et al., (1996) Ester, M., Kriegel, H.-P., Sander, J., and Xu, X. (1996).  A Density-Based Algorithm for Discovering Clusters in Large Spatial Databases with Noise.  Proceedings of the 2nd International Conference on Knowledge Discovery and Data Mining, pages 226–231.
  * Fiksel et al., (2021) Fiksel, J., Datta, A., Amouzou, A., and Zeger, S. (2021).  Generalized bayes quantification learning under dataset shift.  Journal of the American Statistical Association, pages 1–19.
  * Finley et al., (2019) Finley, A. O., Datta, A., Cook, B. D., Morton, D. C., Andersen, H. E., and Banerjee, S. (2019).  Efficient algorithms for bayesian nearest neighbor gaussian processes.  Journal of Computational and Graphical Statistics, 28(2):401–414.
  * Fischer and Modersitzki, (2003) Fischer, B. and Modersitzki, J. (2003).  Curvature based image registration.  Journal of Mathematical Imaging and Vision, 18(1):81–85.
  * Fischer and Modersitzki, (2008) Fischer, B. and Modersitzki, J. (2008).  Ill-posed medicine - An introduction to image registration.  Inverse Problems, 24(3):1–16.
  * Gelman and Rubin, (1992) Gelman, A. and Rubin, D. B. (1992).  Inference from Iterative Simulation Using Multiple Sequences Linked references are available on JSTOR for this article : Inference from Iterative Simulation Using Multiple Sequences.  Statistical Science, 7(4):457–472.
  * Gruen and Baltsavias, (1987) Gruen, A. W. and Baltsavias, E. P. (1987).  High-precision image matching for digital terrain model generation.  Photogrammetria, 42(3):97–112.
  * Grünwald and Mehta, (2020) Grünwald, P. D. and Mehta, N. A. (2020).  Fast rates for general unbounded loss functions: From erm to generalized bayes.  J. Mach. Learn. Res., 21:56–1.
  * Guntupalli et al., (2016) Guntupalli, J. S., Hanke, M., Halchenko, Y. O., Connolly, A. C., Ramadge, P. J., and Haxby, J. V. (2016).  A Model of Representational Spaces in Human Cortex.  Cerebral Cortex, 26(6):2919–2934.
  * Gustin et al., (2012) Gustin, S. M., Peck, C. C., Cheney, L. B., Macey, P. M., Murray, G. M., and Henderson, L. A. (2012).  Pain and plasticity: is chronic pain always associated with somatosensory cortex activity and reorganization?  Journal of Neuroscience, 32(43):14874–14884.
  * Hasson et al., (2004) Hasson, U., Nir, Y., Levy, I., Fuhrmann, G., and Malach, R. (2004).  Intersubject Synchronization of Cortical Activity during Natural Vision.  Science, 303(5664):1634–1640.
  * Haxby et al., (2011) Haxby, J. V., Guntupalli, J. S., Connolly, A. C., Halchenko, Y. O., Conroy, B. R., Gobbini, M. I., Hanke, M., and Ramadge, P. J. (2011).  A common, high-dimensional model of the representational space in human ventral temporal cortex.  Neuron, 72(2):404–416.
  * Hoffman and Gelman, (2014) Hoffman, M. D. and Gelman, A. (2014).  The no-U-turn sampler: Adaptively setting path lengths in Hamiltonian Monte Carlo.  Journal of Machine Learning Research, 15:1593–1623.
  * Iordan et al., (2016) Iordan, M. C., Joulin, A., Beck, D. M., and Fei-Fei, L. (2016).  Locally-optimized inter-subject alignment of functional cortical regions.  arXiv preprint arXiv:1606.02349.
  * Johnson and Christensen, (2002) Johnson, H. and Christensen, G. (2002).  Consistent landmark and intensity-based image registration.  IEEE Transactions on Medical Imaging, 21(5):450–461.
  * Lindquist et al., (2008) Lindquist, M. A. et al. (2008).  The statistical analysis of fmri data.  Statistical science, 23(4):439–464.
  * Lindquist et al., (2009) Lindquist, M. A., Loh, J. M., Atlas, L. Y., and Wager, T. D. (2009).  Modeling the hemodynamic response function in fmri: efficiency, bias and mis-modeling.  Neuroimage, 45(1):S187–S198.
  * Lorbert and Ramadge, (2012) Lorbert, A. and Ramadge, P. J. (2012).  Kernel hyperalignment.  Advances in Neural Information Processing Systems, 3(1):1790–1798.
  * Matheron, (1963) Matheron, G. (1963).  Principles of geostatistics.  Economic geology, 58(8):1246–1266.
  * McAllester, (1999) McAllester, D. A. (1999).  Some pac-bayesian theorems.  Machine Learning, 37(3):355–363.
  * McCarthy et al., (1999) McCarthy, G., Puce, A., Belger, A., and Allison, T. (1999).  Electrophysiological studies of human face perception. ii: Response properties of face-specific potentials generated in occipitotemporal cortex.  Cerebral cortex, 9(5):431–444.
  * Mian et al., (2006) Mian, A. S., Bennamoun, M., and Owens, R. (2006).  Three-dimensional model-based object recognition and segmentation in cluttered scenes.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 28(10):1584–1601.
  * Nenning et al., (2017) Nenning, K. H., Liu, H., Ghosh, S. S., Sabuncu, M. R., Schwartz, E., and Langs, G. (2017).  Diffeomorphic functional brain surface alignment: Functional demons.  NeuroImage, 156(December 2016):456–465.
  * Ombao et al., (2016) Ombao, H., Lindquist, M., Thompson, W., and Aston, J. (2016).  Handbook of Neuroimaging Data Analysis.  CRC Press.
  * Rademacher et al., (1993) Rademacher, J., Caviness Jr, V., Steinmetz, H., and Galaburda, A. (1993).  Topographical variation of the human primary cortices: implications for neuroimaging, brain mapping, and neurobiology.  Cerebral Cortex, 3(4):313–329.
  * Rigon et al., (2020) Rigon, T., Herring, A. H., and Dunson, D. B. (2020).  A generalized bayes framework for probabilistic clustering.  arXiv preprint arXiv:2006.05451.
  * Sabuncu et al., (2010) Sabuncu, M. R., Singer, B. D., Conroy, B., Bryan, R. E., Ramadge, P. J., and Haxby, J. V. (2010).  Function-based intersubject alignment of human cortical anatomy.  Cerebral Cortex, 20(1):130–140.
  * Shawe-Taylor and Williamson, (1997) Shawe-Taylor, J. and Williamson, R. C. (1997).  A pac analysis of a bayesian estimator.  In Proceedings of the tenth annual conference on Computational learning theory, pages 2–9.
  * Simard et al., (2003) Simard, P., Steinkraus, D., and Platt, J. (2003).  Best practices for convolutional neural networks applied to visual document analysis.  In Seventh International Conference on Document Analysis and Recognition, 2003. Proceedings., pages 958–963.
  * Simpson et al., (2012) Simpson, I. J., Schnabel, J. A., Groves, A. R., Andersson, J. L., and Woolrich, M. W. (2012).  Probabilistic inference of regularisation in non-rigid registration.  NeuroImage, 59(3):2438–2451.
  * Thevenaz and Unser, (1998) Thevenaz, P. and Unser, M. (1998).  Efficient mutual information optimizer for multiresolution image registration.  IEEE International Conference on Image Processing, 1:833–837.
  * Thompson et al., (1996) Thompson, P. M., Schwartz, C., Lin, R. T., Khan, A. A., and Toga, A. W. (1996).  Three-dimensional statistical analysis of sulcal variability in the human brain.  Journal of Neuroscience, 16(13):4261–4274.
  * Vehtari et al., (2017) Vehtari, A., Gelman, A., and Gabry, J. (2017).  Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC.  Statistics and Computing, 27(5):1413–1432.
  * Vercauteren et al., (2009) Vercauteren, T., Pennec, X., Perchant, A., and Ayache, N. (2009).  Diffeomorphic demons: Efficient non-parametric image registration.  NeuroImage, 45(1, Supplement 1):S61–S72.  Mathematics in Brain Imaging.
  * Vierck et al., (2013) Vierck, C. J., Whitsel, B. L., Favorov, O. V., Brown, A. W., and Tommerdahl, M. (2013).  Role of primary somatosensory cortex in the coding of pain.  PAIN®, 154(3):334–344.
  * Viola and Wells, (1997) Viola, P. and Wells, W. M. (1997).  Alignment by Maximization of Mutual Information.  International Journal of Computer Vision, 24(2):137–154.
  * Vogt et al., (1995) Vogt, B. A., Nimchinsky, E. A., Vogt, L. J., and Hof, P. R. (1995).  Human cingulate cortex: surface features, flat maps, and cytoarchitecture.  Journal of Comparative Neurology, 359(3):490–506.
  * Vovk, (1990) Vovk, V. G. (1990).  Aggregating strategies.  Proc. of Computational Learning Theory, 1990.
  * Walker and Hjort, (2001) Walker, S. and Hjort, N. L. (2001).  On bayesian consistency.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 63(4):811–821.
  * Watanabe, (2010) Watanabe, S. (2010).  Asymptotic equivalence of Bayes cross validation and widely applicable information criterion in singular learning theory.  Journal of Machine Learning Research, 11:3571–3594.
  * Woo et al., (2015) Woo, C.-W., Roy, M., Buhle, J. T., and Wager, T. D. (2015).  Distinct brain systems mediate the effects of nociceptive input and self-regulation on pain.  PLoS Biol, 13(1):e1002036.
  * Xu et al., (2012) Xu, H., Lorbert, A., Ramadge, P. J., Guntupalli, J. S., and Haxby, J. V. (2012).  Regularized hyperalignment of multi-set fMRI data.  2012 IEEE Statistical Signal Processing Workshop, SSP 2012, pages 229–232.
  * Yang and Medioni, (1992) Yang, C. and Medioni, G. (1992).  Object modelling by registration of multiple range images.  Image and Vision Computing, 10(3):145–155.
  * Yarkoni et al., (2011) Yarkoni, T., Poldrack, R. A., Nichols, T. E., Van Essen, D. C., and Wager, T. D. (2011).  Large-scale automated synthesis of human functional neuroimaging data.  Nature methods, 8(8):665–670.

1 Input: reference image $\mathbf{R}$, reference landmarks $\mathbf{p}^{R}$,
floating image $\mathbf{Y}$, floating landmarks $\mathbf{p}^{Y}$, bound for
distance between landmarks $d$, bound for tuning parameter $\alpha^{max}$.

2 Construct the set
$\mathcal{S}^{p}=\\{\hat{\mathbf{p}}^{Y}\in\mathbf{p}^{Y}:|\hat{\mathbf{p}}^{Y}|=|\mathbf{p}^{R}|\\}$,

3 Construct the corresponding transformation set
$\mathcal{S}^{T}=\\{\hat{T}:\hat{T}=\text{argmin}(\|\mathbf{p}^{R}-T(\hat{\mathbf{p}}^{Y})\|),\text{for}\;\;\hat{\mathbf{p}}^{Y}\in\mathcal{S}^{p}\\}$.

4 for  _each $\alpha\in(0,\alpha^{max})$_ do

5 Find the subset $\hat{\mathcal{S}}^{p}_{\alpha}\subseteq\mathcal{S}^{p}$,
such that
$S(\hat{\mathbf{p}}^{Y}-T_{\alpha}(\hat{\mathbf{p}}^{Y}))+S(\hat{\mathbf{p}}^{R}-T^{-1}_{\alpha}(\hat{\mathbf{p}}^{R}))<\alpha\;\;\text{and}\;\;||\hat{\mathbf{p}}^{R}-T_{\alpha}(\hat{\mathbf{p}}^{Y})||^{2}+||\hat{\mathbf{p}}^{Y}-T^{-1}_{\alpha}(\hat{\mathbf{p}}^{R})||^{2}\leq
2d\\}$

6 Collect the corresponding transformation subset
$\hat{\mathcal{S}}^{T}_{\alpha}\subseteq\mathcal{S}^{T}$,

7 Compute the metric of $\alpha$, $\mathcal{C}(\alpha)$ defined in eqn. (16).

Output: choose $\alpha$ that minimizes $\mathcal{C}(\alpha)$. The
corresponding $T_{\alpha}$ is the proposed prior.

Algorithm 1 Corresponding Feature-based Prior (hyper-parameter) Estimation
Figure 1:  (A) An example of a functional activation map. When used as a
reference map, a query map indicating landmarks of interest, is manually
determined and shown within the dashed rectangle. (B) The query map extracted
from the reference map. (C) The target map generated by warping the reference
map.  Figure 2:  (A) The cropping box comparing the prior information (dashed
box) with the posterior (solid box). (B) A plot of -2 $\times$ expected log
predictive density (ELPD) versus the regularization parameter $\lambda_{T}$
given a fixed $\lambda_{b}$. The vertical bar (red) shows the selected
$\lambda_{T}$, which reaches the minimum of -2$\times$ELPD. Figure 3: The
estimated posterior density plot for each of the transformation parameters:
translations (first row), scalings (second row), and rotations (third row).
Vertical bars indicate the true parameter values. Figure 4: The posterior
means plotted against the true parameters for each transformation parameter:
translations (first row), scalings (second row), and rotations (third row).
Figure 5: Subject-specific activation maps of the region of interest across
the $33$ subjects. Note the differences in location of peaks across subjects,
indicating significant inter-subject differences.  Figure 6: The corresponding
cropped area (bounding box) based on prior information in each subject-
specific map. The features (local peaks) are circled.  Figure 7: The warped
bounding box using the transformation computed with the posterior mean is
shown in solid red. The 95% credible interval of the transformation map is
illustrated in the shaded area.  Figure 8: The warped images obtained using
the diffeomorphic demons registration algorithm implemented in MATLAB’s image
processing toolbox. Note that multiple maps were not registered successfully,
including subjects 3, 5, 6, 12, 19, 25, and 27. Figure 9: A comparison between
the standard approach (left column), the diffeomorphic demons registration
algorithm (middle column), and our proposed approach (right column). Group-
level statistics (t-statistics, group means, group standard deviations) are
calculated for all three approaches. Our proposed approach shows more focal
results and higher t-statistic values, indicating increased sensitivity for
group-level inference.

Figure S1: The posterior means plotted against the true parameters for each
transformation parameter: translations (first row), scalings (second row), and
rotations (third row).