Computational anatomy (CA) is a discipline within medical imaging that studies anatomical shape and form at the visible or gross anatomical scale.
Its scope is broad and draws on anatomy, applied and pure mathematics, medical imaging, neuroscience, physics, probability and statistics, and its focus is the imaged anatomical structure rather than the imaging device.
Within medical imaging, the core concern of computational anatomy is the mapping of information onto anatomical coordinate systems, most often using the dense information in magnetic resonance images (MRI). The introduction of flows into CA parallels the equations of motion in fluid dynamics, and it relies on the idea that dense coordinates in image analysis follow Lagrangian and Eulerian equations of motion. In models based on Lagrangian and Eulerian diffeomorphic flows, the constraints are topological: open sets remain open, coordinates do not cross, which implies uniqueness and the existence of an inverse, and connected sets remain connected. Diffeomorphic methods expanded rapidly and came to dominate the mapping field after the original paper of Christensen and colleagues,[1] and fast and symmetric methods are now available.[2][3]
The core statistical model of computational anatomy in medical imaging is Shannon's source-channel model. The source is a deformable image template I \in \mathcal{I}, and the channel output is the imaging sensor, whose observable is I^D \in \mathcal{I}^{\mathcal{D}}. The importance of the source-channel model is that variation in anatomical configuration and variation in the imaging sensor are modelled separately. Bayesian theory specifies that the model is characterised by a prior \pi_{\mathcal{I}}(\cdot) on the source, acting on I \in \mathcal{I}, together with the conditional density
given I \in \mathcal{I}.
In deformable template theory, an image is related to a template, and the deformation is a group acting on the template. For the image action I(g) \doteq g \cdot I_\mathrm{temp}, g \in \mathcal{G}, a prior \pi_{\mathcal{G}}(\cdot) on the group induces a prior \pi_{\mathcal{I}}(\cdot) on images, and in density form the log posterior is
The random orbit model then specifies how group elements are generated, and therefore how the random objects that make up the prior distribution are drawn.
The random orbit model of computational anatomy first appeared in[4][5][6] and models the coordinate changes produced by group randomness acting on a template, which in turn induces randomness in the image source over orbits of anatomical shape and form and in the observations made by the imaging device. Such a random orbit model, in which randomness on the group induces randomness on images, was studied for object recognition with the special Euclidean group, where the group element g \in \mathcal{G} is an element of that group.[7]
For the study of deformable shapes in CA, the high-dimensional diffeomorphism groups used in computational anatomy are generated by smooth flows \varphi_t, t \in [0,1] that satisfy the Lagrangian and Eulerian flow specifications and the ordinary differential equation
where v \doteq (v_1,v_2,v_3) is a vector field on \mathbb{R}^3 called the Eulerian velocity of a particle of the flow at position \varphi. The vector field is a function in a function space modelled as a smooth Hilbert space, and it is required to have one continuous derivative. For v_t = \dot \varphi_t \circ \varphi_t^{-1}, t \in [0,1], the inverse of the flow is given by
and the 3 \times 3 Jacobian matrix in \mathbb{R}^3 is D\varphi \doteq \left(\frac{\partial \varphi_i}{\partial x_j}\right).
To ensure a smooth flow of diffeomorphisms with inverses, the vector fields on \mathbb{R}^3 must be at least once continuously differentiable in space.[8][9] These vector fields are modelled as elements of a Hilbert space (V, \| \cdot \|_V), and the Sobolev embedding theorem makes each element v_i \in H_0^3, i=1,2,3, have three square-integrable derivatives, so that (V, \| \cdot \|_V) embeds smoothly in the space of continuously differentiable functions.[8][9] The diffeomorphism group consists of flows whose vector fields are absolutely integrable in the Sobolev norm
where \| v_t \|_V^2 \doteq \int_X Av_t \cdot v_t dx, and A is a linear operator A:V \mapsto V^* that defines the RKHS norm. When Av is a generalised function in the dual space V^*, the integral is computed by integration by parts.
In the random orbit model of computational anatomy, the entire flow is reduced to an initial condition that encodes the diffeomorphism. Starting from the initial condition v_0, the geodesic that solves the Euler-Lagrange equation for the computational anatomy Riemannian metric is computed.
Solving the geodesic from the initial condition v_0 is called the Riemannian exponential map, the map from the identity to the group \operatorname{Exp}_\mathrm{id}(\cdot): V \to \operatorname{Diff}_V.
The Riemannian exponential map satisfies \operatorname{Exp}_\mathrm{id}(v_0)= \varphi_1, with initial condition \dot \varphi_0 = v_0 and vector field dynamics \dot \varphi_t = v_t \circ \varphi_t, t \in [0,1].
The map extends to the whole group as \varphi= \operatorname{Exp}_\varphi(v_0\circ \varphi) \doteq \operatorname{Exp}_\mathrm{id}(v_0) \circ \varphi.
A random orbit is generated by randomising the flow, that is, by drawing an initial tangent vector field v_0 \in V at the identity, to produce a random object n \doteq \operatorname{Exp}_\mathrm{id}(v_0) \cdot m_0 \in \mathcal{M}.
Conditioned on a particular atlas I_a \in \mathcal{I}, the random orbit model induces a prior on shapes and images I \in \mathcal{I}. The generative model produces a mean field I as a random change of template coordinates, I \doteq \varphi \cdot I_a, where the diffeomorphic coordinate change is generated randomly by a geodesic flow.
Conditioned on a particular atlas I_a \in \mathcal{I}, the random orbit model induces a prior on shapes and images I \in \mathcal{I}. The generative model produces a mean field I as a random change of template coordinates, I \doteq \varphi \cdot I_a, with the diffeomorphic coordinate change generated by a geodesic flow. The prior \pi_\mathrm{Diff}(d\varphi) on \operatorname{Diff}_V is induced by the flow \operatorname{Exp}_\mathrm{id}(v), where v \in V is constructed from a Gaussian random field prior \pi_V(dv). The density of the random observable at the sensor output, I^D \in \mathcal{I}^D, is
Maximum a posteriori (MAP) estimation is central to modern statistical theory. The parameter of interest \theta \in \Theta can take several forms, including (i) disease type, such as neurodegenerative or neurodevelopmental disease, (ii) structure type, such as cortical or subcortical structures relevant to image segmentation, and (iii) reconstruction of a template from a population. Given an observed image I^D, the MAP estimate maximises the posterior
This requires the conditional probability p(\theta\mid I^D) = \frac{p(I^D,\theta)}{p(I^D)}. The multi-atlas orbit model randomises over a countable set of atlases \{ I_a, a \in \mathcal{A} \}. The model of images in the orbit takes the form of a multimodal mixture
Conditional Gaussian models have been studied extensively for dense image matching and for landmark matching.
The image I^D(x), x \in X is modelled as a conditional Gaussian random field with conditional mean field \varphi_1 \cdot I \doteq I(\varphi_1^{-1}), \varphi_1 \in \operatorname{Diff}_V. For uniform variance, the endpoint error term plays the role of the log conditional probability, depending only on the mean field, and gives the endpoint term.
The landmarks Y = \{ y_1,y_2,\dots \} are modelled as conditional Gaussian with mean field \varphi_1(x_i), i=1,2,\dots, \varphi_1 \in \operatorname{Diff}_V, and with noise variance independent of the landmark and constant. The log conditional probability, which depends only on the mean field, is treated as the endpoint term:
The random orbit model of multiple atlases models the orbit of a shape as the union of several anatomical orbits generated by the action of the diffeomorphism group, \mathcal{I} = \bigcup_{a \in \mathcal{A}}\operatorname{Diff}_V \cdot I_a. Each atlas has a template and a predefined segmentation field (I_a,W_a), a = a_1,a_2,\ldots that partitions the anatomical structures in MRI coordinates. These pairs are indexed on a voxel grid I_a(x_i),W_a(x_i), x_i \in X \subset \mathbb{R}^3, and they contain an MRI image and a dense label at each voxel coordinate. The anatomical labels of the partitioned structures are drawn manually by neuroanatomists.
The Bayesian segmentation problem[10] given measurements I^D, with mean field and partition (I,W), requires estimation of the anatomical labels \theta \doteq W. The mean field of the observable image I^D is modelled as a random deformation of some template I \doteq \varphi \cdot I_a, and the template is itself randomly selected, A = a. The optimal diffeomorphism \varphi \in \mathcal{G} is hidden, and it acts on the coordinate background space of the randomly selected template image I_a. Given a single atlas a, the likelihood model for inference is determined by the joint probability p(I^D,W\mid A = a). For multiple atlases, fusion of the likelihood functions produces a multimodal mixture model, and the prior averages over models.
The MAP estimator of the segmentation W_a maximises \max_{W} \log p(W \mid I^D) with I^D given, over the mixture of all atlases:
The quantity p(I^D,W) is computed by likelihood fusion across deformable atlases, where \pi_A(a) is the prior probability that the observed image comes from a particular template image I_a.
MAP segmentation can be solved iteratively with the expectation-maximisation algorithm:
Generating a template empirically from a population is a basic operation throughout the discipline.
Several Bayesian methods have appeared for submanifolds and for dense image volumes.
For the dense image volume case, given observations I^{D_1}, I^{D_2}, \dots, the problem is to estimate a template in the dense image orbit I \in \mathcal{I}. Ma's procedure starts from an initial hypertemplate I_0 \in \mathcal{I} and models the template in the orbit through an unknown diffeomorphism to be estimated, I \doteq \varphi_0 \cdot I_0, with parameter the log coordinate \theta \doteq v_0 that determines the geodesic mapping of the hypertemplate \operatorname{Exp}_\mathrm{id}(v_0) \cdot I_0 = I \in \mathcal{I}.
In the Bayesian random orbit model of computational anatomy, the observed MRI images I^{D_i} are modelled as conditional Gaussian random fields with mean field \varphi_i \cdot I, where \varphi_i is a random unknown transformation of the template. The MAP estimation problem estimates the unknown template I \in \mathcal{I} from the observed MRI images.
Ma's dense image procedure starts from an initial hypertemplate I_0 \in \mathcal{I} and models the template in the orbit through an unknown diffeomorphism to be estimated, I \doteq \varphi_0 \cdot I_0. The observations are modelled as conditional random fields, with I^{D_i} a random field with mean field \varphi_i \cdot I \doteq \varphi_i \cdot \varphi_0 \cdot I_0. The unknown variable explicitly estimated by MAP is the hypertemplate mapping \varphi_0, and the other mappings are treated as nuisance or hidden variables and integrated out by the Bayesian procedure. This is done with the expectation-maximisation algorithm.
The orbit model is exploited by associating the unknown flows to be estimated with their log coordinates v_i,i=1,\dots, using the logarithm and exponential maps of Riemannian geodesics in computational anatomy, that is, the initial vector fields in the tangent space at the identity, so that \operatorname{Exp}_\mathrm{id}(v_i) \doteq \varphi_i, where \operatorname{Exp}_\mathrm{id}(v_0) is the mapping of the hypertemplate.
The MAP estimation problem becomes
The EM algorithm treats the complete data as the vector field coordinates of the parameterised mappings v_i,i=1,\dots, and iteratively computes the conditional expectation
References