An Algorithm for Learning Shape and Appearance Models without Annotations
John Ashburner Mikael Brudfors Kevin Bronik Yael Balbastre



mymainaddress]John Ashburner
j.ashburner@ucl.ac.uk
mymainaddress]Mikael Brudfors
mymainaddress]Kevin Bronik
mymainaddress]Yaël Balbastre


[mymainaddress]
Wellcome Centre for Human Neuroimaging

UCL Queen Square Institute of Neurology

12 Queen Square, London, WC1N 3BG, UK.



This paper presents a framework for automatically learning shape and appearance models for medical (and certain other) images.
It is based on the idea that having a more accurate shape and appearance model leads to more accurate image registration, which in turn leads to a more accurate shape and appearance model.
This leads naturally to an iterative scheme, which is based on a probabilistic generative model that is fit using Gauss-Newton updates within an EM-like framework.
It was developed with the aim of enabling distributed privacy-preserving analysis of brain image data, such that shared information (shape and appearance basis functions) may be passed across sites, whereas latent variables that encode individual images remain secure within each site.
These latent variables are proposed as features for privacy-preserving data mining applications.

The approach is demonstrated qualitatively on the KDEF dataset of 2D face images, showing that it can align images that traditionally require shape and appearance models trained using manually annotated data (manually defined landmarks etc.).
It is applied to MNIST dataset of handwritten digits to show its potential for machine learning applications, particularly when training data is limited.
The model is able to handle “missing data”, which allows it to be cross-validated according to how well it can predict left-out voxels.
The suitability of the derived features for classifying individuals into patient groups was assessed by applying it to a dataset of over 1,900 segmented T1-weighted MR images, which included images from the COBRE and ABIDE datasets.



Machine Learning Latent Variables Diffeomorphisms Geodesic Shooting Shape Model Appearance Model











§ INTRODUCTION

This paper introduces an algorithm for learning a model of shape and appearance variability from a collection of images, without relying on manual annotations.
The shape part of the model concerns modelling variability with diffeomorphic deformations, which is essentially image registration.
In contrast, the appearance part is about accounting for signal variability that is not well described by deformations, and is essentially about adapting a “template” to enable more precise registration.

The problem of image registration is often viewed from a Bayesian perspective, whereby the aim is to determine the most probable deformation (ψ) given the fixed (f) and moving (μ) images 

    ψ̂=_ψlog p(ψ |  f, μ) = _ψ( log p( f | ψ, μ) + log p(ψ) ).

In practice, the regularisation term (log p(ψ)) is not usually defined empirically, and simply involves a penalty based on some simple measure of deformation smoothness.
One of the aims of this work is to try to improve on this simple model.
By providing empirically derived priors for the allowable deformations, trained shape models have been shown to exhibit more robust image registration.
An early example is <cit.>, in which control point positions are constrained by their first few modes of variability.
Training this model involved annotating images by manually placing a number of corresponding landmarks, computing the mean and covariance of the collection of landmarks, and then computing the eigenvectors of the covariance <cit.>.
In neuroimaging, shape models have previously been used to increase the robustness of brain image segmentation <cit.>.
The current work involves densely parameterised shape models within the diffeomorphic setting, and relates to previous work on diffeomorphic shape models <cit.>, as well as those using more densely parameterised deformations <cit.>.
Recently, <cit.> developed their Principal Geodesic Analysis (PGA) framework for directly computing the main modes of shape variation within a diffeomorphic setting.
The work presented in this paper borrows heavily from that of <cit.>, but extends it to involve a joint model of shape and appearance variability.

In addition to increasing the robustness of image registration tasks, shape models can also provide features that may be used for statistical shape analysis.
This is related to approaches used in geometric morphometrics <cit.>, where the aim is to understand shape differences among anatomies.
Shape descriptors from the PGA framework have previously been found to be useful features for data mining <cit.>.


A number of previous works have investigated combining both shape and appearance variability into the same model <cit.>.
These combined shape and appearance models have generally shown good performance in a number of medical imaging challenges <cit.>.
While there is quite a lot written about learning appearance variability alone, the literature on automatically learning both shape and appearance together is fairly limited. 
Earlier approaches required annotated data for training, but there are now some works appearing that have looked into the possibility of using unsupervised or semi-supervised approaches for learning shape and appearance variability.
Examples include <cit.>, <cit.>, <cit.> and <cit.>.
The current work is about an unsupervised approach, but there is no reason why it could not be made semi-supervised by also incorporating some manually defined landmarks or other features.

This work was undertaken as a task in the Medical Informatics Platform of the EU Human Brain Project (HBP).
The original aim of the Medical Informatics Platform was to develop a distributed knowledge discovery framework that enables data mining without violating patient confidentiality.
The strategy was to involve a horizontally partitioned dataset, where data about different patients is stored in different hospital sites, and only data aggregated over many subjects may leave any site.
While aggregates could be exploited by those with malicious intent in order to extract some information about individual patients, it is generally more difficult to do so, particularly when there are constraints on how much aggregated data may leave each site.
The algorithm presented in this paper can be implemented in a way that does not require patient-specific information to leave a site, and instead only shares aggregates, which reveal less about the individual subjects.
Some leakage of information is inevitable, particularly for sites holding data on only small numbers of individuals, but we leave this as a topic to be addressed elsewhere.

For data mining, in situations where the dimensionality of the data (e.g. number of voxels in an image) is greater than the number of data points (e.g. number of images), it is sometimes possible to make use of the Woodbury matrix identity[The formula is frequently encountered as (Σ^-1+ F S^-1 F^T)^-1 = Σ - Σ F( S +  F^TΣ F)^-1 F^TΣ (see <https://en.wikipedia.org/wiki/Woodbury_matrix_identity>). Typically, F is an M × N feature matrix, Σ is an M × M prior covariance and S is an N × N matrix (usually diagonal) that encodes the residual covariance. This is useful when M ≫ N, as inverses of only N × N matrices are required. Chapter 2 of <cit.> shows how this relates to kernel methods.], which leads naturally to kernel-based approaches for machine learning.
However, these would require dot-products or differences to be computed between data in different sites, and this would be prohibited by the privacy preserving framework.
Also, the Woodbury identity becomes less useful for extremely large datasets, because of both memory requirements and computational complexity.
Aggregated data may be weighted moments (e.g. ∑_n r_n, ∑_n r_n  z_n or ∑_n r_n  z_n  z_n^T, where z_n is a vector of values for patient n, and r_n is a patient-specific weight generated by some rule), which could then be used for clustering or other forms of statistical analysis.
Dimensionality reduction is often needed to enable effective data mining of images, particularly if covariances need to be represented (such as for clustering into patient subgroups using Gaussian mixture models).

Recently, methods such as convolutional neural networks have begun to show a great deal of promise for machine learning tasks <cit.>, and it seems that much of the previous work from the field of medical imaging will be replaced by new ideas from the machine learning field.
This is partly due to the very large sets of labelled training data now available for certain applications.
When the number of training samples is large compared to the dimensionality of each example, these methods can accurately capture the non-linearities in the training data.
Recent MICCAI and other medical imaging conferences have shown that this is especially true for tasks such as image segmentation, where each voxel is effectively a label.
However, they can be less accurate when the number of training samples is relatively small compared to the complexity of each data point, when other approaches, such as image registration, may be better at encoding the non-linearities <cit.>.
Conventional deep learning simply models the output as a nonlinear function of the input.
In contrast, generative approaches to machine learning would involve constructing probabilistic models of the input, such that random samples could be drawn that are as similar as possible to possible future examples of real input data.
A recent incarnation is the family of generative adversarial networks <cit.>, where the general aim is to learn a model that can generate samples that are indistinguishable from real data.
Accurate generative models of data have many potential applications, which include outlier detection and augmenting training data for deep learning.
There is increasing interest in the use of generative approaches for machine learning, partly because they can be extended to work in a semi-supervised way.
This enables unlabelled training data to contribute towards the model, potentially allowing more complex models to be learned from fewer labelled examples.
A simple example of semi-supervised learning would be the approach of <cit.>, which used both labelled and unlabelled images for generating tissue probability maps.

Another motivation for generative modelling approaches is to enable missing data to be dealt with.
Brain images – particularly hospital brain images – often have different fields of view from each other, with parts of the brain missing from some of the scans.
Many machine learning approaches do not work well in the presence of missing data, so imputing missing information is an implicit part of the framework presented in this paper.


This work proposes a solution based on learning a form of shape and appearance model.
The overall aim is to capture as much anatomical variability as possible using a relatively small number of latent variables.
In addition to 3D brain image data, a number of other types of images will be used to illustrate other aspects of the very general framework that we present.
Often, it can be difficult to understand the behaviour of an algorithm when it is only applied to 3D volumes, particularly when space for figures in papers is limited.
For this reason, the work is also applied to a number of 2D images, ranging from the tiny images that make up the MNIST dataset, to some images of faces, along with some single slices through brain images.







§ METHODS

The proposed framework builds heavily on the principal geodesic analysis model of <cit.>.
Modifications involve extending the framework to use a Gauss-Newton optimisation strategy and incorporating an appearance model.
The Gauss-Newton strategy is intended to be more suited to the distributed computing framework of the HBP, although a streaming method using stochastic gradient descent could also have been used.

The basic idea is that both shape and appearance may be modelled by linear combinations of spatial basis functions, and the objective is to automatically learn the best set of basis functions and latent variables from some collection of images.
This is essentially a form of factorisation of the data.
Each of the N images will be denoted by f_n ∈ℛ^M, where M is the number of pixels/voxels in an image, 1≤ n≤ N, and the entire collection of images by F.
An appearance model for the nth image is constructed from a linear combination of basis functions, such that

    a_n =μ +  W^a  z_n.

Here, W^a is a matrix containing K appearance basis functions, and z_n is a vector of K latent variables for the nth image.
The vector μ is a mean image, with the same dimensions as a column of W^a.
The shape model is encoded similarly, with

    v_n = W^v  z_n.

The approach involves warping a_n to match f_n, using a diffeomorphic deformation parameterised by v_n.
The diffeomorphism (ψ_n) is constructed from v_n by a procedure known as “geodesic shooting”, which is outlined in Section <ref>.
From a probabilistic perspective, the likelihood can be summarised by

    p( f_n |   z_n, μ,  W^a,  W^v) = p( f_n |  a_n(ψ_n)).

Different forms of appearance model are presented in Section <ref>, but for convenience, we use the generic definition

    J( f_n,  z_n, μ,  W^a,  W^v) = -ln p( f_n |   z_n, μ,  W^a,  W^v).


In practice, a small amount of regularisation is imposed on the mean (μ) by assuming it is drawn from a multivariate Gaussian distribution of precision L^μ

    p(μ) =   𝒩(μ |  0, ( L^μ)^-1).

The precision matrix is only involved conceptually.
Because the matrix is circulant (a special kind of Toeplitz matrix where data are assumed to wrap around at the boundaries), matrix multiplication and division may be effected by 3D convolutions.
Possible forms for L^μ are described in Section <ref>.

A weighted average of two strategies for regularising the basis functions (W^a and W^v) and latent variables (z_n) is used, which are:


  * The first strategy involves separate priors on the basis functions, and on the latent variables.
Each of the basis functions is assumed to be drawn from zero-mean highly multivariate Gaussian, parameterised by very large and sparse precision matrices.
Possible forms of the matrices for regularising shape (L^v) are described in Section <ref>, whereas those for appearance (L^a) are described in Section <ref>.
Priors for the basis functions are

    p( W^v) =   ∏_k=1^K^v𝒩( w^v_k |  0, (N  L^v)^-1), 
    
    p( W^a) =   ∏_k=1^K^a𝒩( w^a_k |  0, (N  L^a)^-1).

The latent variables (Z) are assumed to be drawn from zero-mean multivariate Gaussian distributions, parameterised by a precision matrix (A) that is derived from the data[Note that the precision matrix A should not be confused with the variables a_n, which were introduced earlier. Hopefully, the context in which they are used should be enough to prevent any confusion.].

    p( z_n |  A) =𝒩( z_n |  0,  A^-1).

Note that precision matrices for the basis functions, L^v and L^a, are scaled by N.
The reason for this is that it keeps the variance of the latent variables within a similar range, irrespective of how many images are involved, which makes it easier to define priors for them.

The simplest approach for working with the A matrix is to make a point estimate based on the distribution of the latent variables.
This approach allows the relevance of each of the basis functions (columns of W^a and W^v) to be automatically determined.
If they are not needed, some of the elements of the latent variables (i.e. rows of Z) are forced to zero <cit.>.
In practice, this could lead to prematurely sparse solutions that may be difficult to recover from.
Instead, sparsity is avoided by assuming that matrix A is drawn from a Wishart distribution.

    p( A) =   𝒲_K( A | Λ_0, ν_0)=   | A|^(ν_0 - K - 1)/2exp (-12(Λ_0^-1 A))/2^(ν_0 K)/2 |Λ_0|^ν_0/2Γ_K(ν_02),

where Γ_K is the multivariate gamma function.
This prior can be made as uninformative as possible by using ν_0 = K and Λ_0 =  I/ν_0, where I is an identity matrix.
In general, Λ_0 should be a positive definite symmetric matrix, with ν_0 ≥ K so that the distribution can be normalised.


  * The second strategy is similar to that of <cit.>, and is a pragmatic solution to ensuring that enough regularisation is used.

    ln p( Z, W^a, W^v) = -12( Z Z^T (( W^a)^T  L^a  W^a + ( W^v)^T  L^v  W^v)) + const

This strategy imposes smoothness on the reconstructions by assuming penalties based on ln𝒩( W^a  z_n |  0,  L^a) and ln𝒩( W^v  z_n |  0,  L^v), in a similar way to more conventional regularisation approaches.


The relative weighting of the two strategies is by λ_1 and λ_2.
Better results are typically achieved when hyper-parameters are specified so that they sum to a value greater than 1.
When everything is combined (see Fig. <ref>), the following joint log-probability is obtained

    ln p( F,    μ,  W^a,  W^v,  A,  Z)=    - ∑_n=1^N J( f_n,  z_n, μ,  W^a,  W^v) - 12μ^T  L^μμ    - λ_1 N2( (( W^a)^T  L^a  W^a) + (( W^v)^T  L^v  W^v) )    + λ_12((N + ν_0 - K - 1) ln | A| - (( Z Z^T + Λ_0^-1)  A) )    - λ_22( Z Z^T (( W^a)^T  L^a  W^a + ( W^v)^T  L^v  W^v)) + const.




Fitting the model is described in Appendix <ref>.
Ideally, the procedure would compute distributions for all variables, such that uncertainty was dealt with optimally.
Unfortunately, this is computationally impractical for the size of the datasets involved.  Instead, only point estimates are made for the latent variables (ẑ_n) and various parameters (μ̂, W^a, W^v, and sometimes σ̂^2), apart from A, which is inferred within a variational Bayesian framework.


The approach also allows shapes and appearances to be modelled separately, by having some of the latent variables control appearance, and others control shape.
This may be denoted by

    a_n =   μ + ∑_k=1^K^a w^a_k z_k,n,
     v_n =   ∑_k=1^K^v w^v_k z_K^a + k, n.

For simplicity, only the form where each latent variable controls both shape and appearance is described in detail.
This is the form used in active appearance models <cit.>.
Note however, that in the form where shape and appearance are controlled by separate latent variables, the precision matrix A still encodes covariance between the two types of variables.  This means that latent variables controlling either shape or appearance are not estimated completely independently.




 §.§ Shape Model 

Relative shapes are encoded by diffeomorphisms, which are each parameterised by a set of initial conditions (v_n =  W^v  z_n).
This subsection mixes both discrete and continuous representations of the same objects, so some notation is now introduced.
For the discrete case, where a velocity field is treated as a vector, it are denoted by v_n.
Alternatively, the same object may be treated as a continuous 3D vector field, where it is denoted by v_n.

In addition, deformations may be treated as discrete or continuous.
Within the continuous setting, warping an image by a diffeomorphism ψ may be denoted by a' = a(ψ).
In the discrete setting, this resampling may be conceptualised as a matrix multiplication, where a very large sparse matrix Ψ encodes the same deformation (and associated trilinear interpolation), such that a' = Ψ a.
The transpose of this matrix can be used to perform a push-forward operation, which is frequently used in this work and which we denote by f' = Ψ^T  f.

Diffeomorphisms are created within the Large-Deformation Diffeomorphic Metric Mapping (LDDMM) framework <cit.>, which allows image registration to produce smooth, invertible one-to-one mappings.
A simplification of the basic idea is that large deformation diffeomorphisms may be computed by composing together a series of much smaller deformations.
Providing the constituent deformations are one-to-one, then the result of their composition should also be one-to-one.
The original LDDMM implementation involved a variational optimisation, which involves estimating a series of velocity fields (v_t) at different “time-points”, t.
An alternative to the variational approach involves estimating only an initial velocity (v_0), and deriving velocities at subsequent time points by a procedure known as “geodesic shooting” <cit.>.
Algorithm <ref> in the Appendix presents a simple algorithm for geodesic shooting, which generates a deformation ψ from an initial velocity field v.

Another important feature is that when image registration is formulated as a generative model of the individual images, the diffeomorphic framework can be set up such that the parametrisation of the deformations (i.e., the initial velocity fields) are all in alignment with each other in the space of the template.
This makes them more useful for various forms of multivariate models of shape variability as they are better encoded by factorisation models.




  §.§.§ Differential operator for shape model 
 
Our implementation of geodesic shooting uses Fast Fourier Transform (FFT) methods to obtain the Green's function <cit.>.
Because FFTs are used throughout the work, the boundary conditions for the velocity fields are assumed to be periodic.
The precision matrix used in Eqn. <ref> has the form

    v^T  L^v  v=   ∫_x∈Ω(ω^v_0  v(x) ^2 + ω^v_1 ∇ v(x)^2 + ω^v_2 ∇^2 v(x)^2 ) dx
     +    ∫_x∈Ω(ω^v_3/4 Dv(x) + (Dv(x))^T_F^2  + ω^v_4 (Dv(x))^2 
       ) dx

where |·|_F denotes the Frobenius norm (the square root of the sum of squares of the matrix elements).
The above integral is defined in Sobolev space, which is a weighted Hilbert space where spatial derivatives, up to a certain degree, are accounted for.

Five hyper-parameters are involved:


  * ω^v_0 controls absolute displacements, and is typically set to be a very small value.

  * ω^v_1 controls stretching, shearing and rotation.

  * ω^v_2 controls bending energy.  This ensures that the resulting velocity fields have smooth spatial derivatives.

  * ω^v_3 controls stretching and shearing (but not rotation).

  * ω^v_4 controls the divergence, which in turn determines the amount of volumetric expansion and contraction.

Most of the regularisation in this work was typically based on a combination of the linear-elasticity (ω^v_3 and ω^v_4) and bending energy (ω^v_2) penalties.
The differential operator involved very little penalty against absolute displacements (ω^v_0), although some was necessary to ensure that the Green's function could be computed, while still allowing the deformations to incorporate uniform translations.
The types of differential operators used in <cit.> would be constructed from ω^v_0, ω^v_1 and ω^v_2, whereas those used in <cit.> were largely based on Lamé's constants ω^v_3 and ω^v_4.
The latter parametrisation allows the stresses and strains of deformations to be formulated in terms of, for example, Poisson's ratio and Young's modulus[See <https://en.wikipedia.org/wiki/Linear_elasticity>.].




 §.§ Appearance Models 

A number of different choices for the appearance model are available for Eqn. <ref>, each suitable for modelling different types of image data.
These models are based on p( f_n |  a'_n), which leads to an “energy” term (J) that drives the model fitting and is assumed to be independent across voxels

    a'_n    =Ψ_n (μ +  W^a  z_n)
    
    J( a'_n)    = -ln p( f_n |   a'_n) = -∑_m=1^M ln p(f_mn | a'_mn).

Because the approach is generative, missing data can be handled by simply ignoring those voxels where there is no information.
By doing this, they do not contribute towards the objective function and play no role in driving the model fitting.
A number of different “energy” functions have been implemented for modelling different types of data.
These are listed next.




  §.§.§ Gaussian noise model

Mean-squares difference is a widely used objective functions for image matching, which is based on the assumption of stationary Gaussian noise.
For an image consisting of M pixels or voxels, the function would be

    -J_L_2( a') =ln p( f |  a', σ^2) = -M2ln(2π) - M2lnσ^2 - 12σ^2|| f- a'||_2^2,

where ||·||_2 denotes the Euclidean norm.
The simplest approach to compute σ^2 is to make a maximum likelihood estimate from the variance by

    σ̂^̂2̂=1MN∑_n=1^N || f_n- a'_n||_2^2 .





  §.§.§ Logistic function with Bernoulli noise model

When working with binary images, such as single tissue type maps having voxels of zeros and ones (or values very close to zero or one), it may be better to work under the assumption that voxels are drawn from a Bernoulli distribution, which is a special case of the binomial distribution.
For a single voxel,

    P(f | s) = s^f(1-s)^1-f.


The range 0<s<1 must be satisfied, which can be achieved using a logistic sigmoid function

    s(a') =1/1+exp(-a').


Putting these together gives

    P(f | a')    =(1/1+exp(-a'))^f(1-1/1+exp(-a'))^1-f   =exp(a' f) s(-a').


This gives the matching function

    -J_Bern( a') =ln P( f |  a') = ∑_m=1^M (f_m a'_m + ln s(-a'_m)).





  §.§.§ Softmax function with categorical noise model

If there are several binary maps to align simultaneously, for example maps of grey matter, white matter and background, then a categorical noise model is appropriate.
A categorical distribution is a generalisation of the Bernoulli distribution, and also a special case of the multinomial distribution.
The probability of a vector f of length C, such that f_c ∈{0,1} and ∑_c=1^C f_c = 1, is given by 

    P( f |  s) =∏_c=1^C s_c^f_c,

where s_c>0 and ∑_c=1^C s_c = 1.
The constraints on s can be enforced by using a softmax function.

    s_c( a') =exp a'_c/∑_c=1^C exp a'_c

Using the “log-sum-exp trick”, numerical overflow or underflow can be prevented by first subtracting the maximum of a, so

    s_c( a') =exp(a'_c-a^*)/∑_c=1^C exp(a'_c-a^*), where  a^* =max{a'_1,,a'_C}


Noting that each image is now a matrix of M voxels and C classes, the objective function can then be computed as

    -J_cat( A')    =ln P( F |  A')   =∑_m=1^M (∑_c=1^C a'_mc f_mc - a^* - log(∑_c=1^C exp(a'_mc-a_m^*) ) )





  §.§.§ Differential operator for appearance model 

Regularisation is required for the appearance variability, as it helps to prevent the appearance model from absorbing too much of the variance, at the expense of the shape model.
This differential operator (again based on a Sobolev space) is used in Eqns. <ref> and <ref>, and controlled by three hyper-parameters.

    a^T  L^a  a=∫_x ∈Ω(ω^a_0  a(x) ^2 + ω^a_1 ∇ a(x)^2 + ω^a_2 ∇^2 a(x)^2 ) dx

The extent to which intensity differences may influence apparent spatial deformations was briefly discussed in <cit.>.
Automatically disambiguating between what should be modelled as a shape change, and what should be treated as as appearance change, is still a largely unresolved area.
The approach taken in this work is simply to treat the construct as a model of the data, and to assess it according to how well it describes and predicts the observations, rather than how well it can separately estimate shape information versus appearance information.







§ RESULTS

To show the general applicability of the approach, evaluations were performed with a number of datasets of varying characteristics.
Our implementation[Available from <https://github.com/WTCN-computational-anatomy-group/Shape-Appearance-Model>.] is written in a mixture of MATLAB and C code (MATLAB “mex” files for the computationally expensive parts).




 §.§ Qualitative 2D experiments with faces

After years of exposure to faces, most people can identify whether an image of a face is plausible or not, so images of human faces provide a good qualitative test of how well the algorithm can model biological variability.
Recent evidence <cit.> suggests that the primate brain might even use some form of shape and appearance model for encoding faces, which might indicate that the type of model used in this work could be quite effective.

The straight on views from the Karolinska Directed Emotional Faces (KDEF) data-set <cit.> were used to make a visual assessment of how well the algorithm performs.
This data-set consisted of photographs of 70 participants, holding seven different facial expressions, which was repeated twice.  Some of the images were not usable, so the final dataset consisted of 932 colour images, which were downsampled to a size of 282 × 382.
The original intensities were in the range of 0 to 255, but these values were re-scaled by 1/255.

A 64 eigenmode model was used (K=64), which assumed Gaussian noise.
Model fitting (i.e., learning the shape and appearance basis functions, etc.) was run for 20 iterations, with ν_0 = 1000, λ = [15.2  0.8], ω^a = [4  512  64], ω^μ = N [10^-4 0.1  0.1] and ω^v = [10^-3 0  16  1  1].
It was fit to the entire field of view of the images, rather than focusing only on the faces, and some of the resulting fits are shown in Fig. <ref>.
The first set of images are a random selection of the original data, with the full shape and appearance model fits shown immediately below.
As can be seen, the fit is reasonably good - especially given that only 64 modes of variability were used, and that these have to account for a lot of variability of hair etc.
Below these are the shape model fits, generated by warping the mean according to the estimated deformations (μ(ψ_n)).
The appearance fits are shown at the bottom (a_n from Eqn. <ref>).
Ideally, these reconstructions of appearance should be in perfect alignment with each other, which is not quite achieved in certain parts of the images.
In particular, the thickness of the neck varies according to whether or not the people in the images have short or long hair.
When looked at separately, the shape and appearance parts of the model do not behave quite so well, but when combined, they give quite a good fit.
Fig. <ref> shows a simple 64-mode principal component analysis (PCA) fit to the same data, which clearly does not capture variability quite as well as the shape and appearance model.





The first 12 modes of variability are illustrated in Fig. <ref>.
Much of the variance to be modelled is not actually part of the face, which can be seen in how the hairstyles of the subjects drive the first few modes.
Better fits would have been achieved if only the faces themselves had been modelled, but the main aim here was to illustrate the behaviour of the approach in more challenging situations.
Masking out parts of the data from model fitting is possible, as one of the objectives of the work was to be able to handle missing data (see Section <ref>).




For these examples, there should really have been a distinction between inter-subject variability and intra-subject variability, using some form of hierarchical model for the latent variables.
This type of hierarchical mixed-effects model is widely used for analysing multi-subject data within the neuroimaging field <cit.>, and a number of works have applied mixed effects modeling to image registration <cit.>.

Spatial correlations in the images mean the assumptions about the independence of the noise across neighbouring voxels to be less valid, which causes problems for the simple i.i.d. noise model used here.
This is why the regularisation needed to be weighted more heavily relative to the image matching term.
This is related to the “virtual decimation” <cit.> approach of down-weighting the matching term by a correction factor that accounts for the number of independent observations.
Numbers of independent observations could be derived using random field theory <cit.>, using gradients of residuals to estimate image smoothness in terms of the full-width half maximum of a Gaussian.




  §.§.§ Simulating faces

Once the model is learned, it becomes possible to generate random faces from the estimated distribution.
This involves drawing a random vector of latent variables z∼𝒩( 0,Â^-1), and using these to reconstruct a face.
Fig. <ref> shows two sets of randomly generated faces, where the lower set used the same latent variables as the upper set, except that they were multiplied by -1.
Although some of the random faces are not entirely plausible, they are much more realistic than faces generated from a simple 64-mode PCA model (shown in Fig. <ref>).






  §.§.§ Vector arithmetic

In many machine learning applications, it is useful to be able to model certain non-linearities in the data in a linear way, allowing more interpretable linear methods to be used while still achieving a good fit.
Following <cit.>, this section shows that simple arithmetic on the latent variables can give intuitive results.
The first three columns of Fig. <ref> show the full shape and appearance model fits to various faces.
Images in the right hand column of Fig. <ref> were generated by making linear combinations of the latent variables that encode the images in the first three columns, and then reconstructing from these.
Unlike arithmetic computed in pixel space (not shown), performing arithmetic on the vectors encoding the images gives reasonably plausible results.
For medical imaging applications, it may be useful (for example) to transport disease-related changes computed from one patient, or set of patients, on to the image of a new individual <cit.>.
Even when not explicitly transporting such information, it is useful to have models where the encoding can be treated in an approximately linear way.


 



 §.§ 2D experiments with MNIST

In this section, the behaviour of the approach using “big data” is assessed, which gives more of an idea of how this type of method may behave with some of the very large image datasets currently being collected.
Instead of testing on a large collection of medical images, the approach was applied to a large set of tiny images of hand-written digits.
MNIST[<http://yann.lecun.com/exdb/mnist/>.] <cit.> is a modified version of the handwritten digits from the National Institute of Standards and Technology (NIST) Special Database 19.
The dataset consists of a training set of 60,000 28 × 28 pixel images of the digits 0 to 9, along with a testing set of 10,000 digits.
MNIST has been widely used for assessing the accuracy of machine learning approaches, and is used here as it allows behaviour of the current approach to be compared against the state-of-the-art pattern recognition methods.

In recent years, the medical imaging community has seen many of the established “old-school” approaches replaced by deep learning, but in doing so, “have we thrown the baby out with the bath water?”[This was said by the late David MacKay <cit.> in relation to the success of kernel methods, such as support-vector machines or Gaussian processes, which, at the time, were replacing neural networks in practical applications.].
There may still be widely used concepts from orthodox medical imaging (i.e., not deep learning) that are still useful.
In particular, geometric transformations of images are now finding their way into various machine learning approaches (e.g. <cit.>).
Much of the early work on deep learning was performed using MNIST.
Although good accuracies were achieved, the computer vision community did not take such work seriously because the images were so small.
This, however, was the early days of deep learning (i.e., before 2012), and was a sign of things to come.
This section describes an attempt to begin to reclaim some of the territory lost to deep learning.

Unlike most conventional pattern recognition approaches, the strategy adopted here is generative.
Training involves learning independent models of the ten different digits in the training set, while testing involves fitting each model in turn to each image in the test set, and performing model comparison to assess which of the ten models better explains the data.
The training stage involved learning μ̂, Ŵ^a, Ŵ^v and Â for each digit class.
A similar strategy was previously adopted by <cit.>.
From a probabilistic perspective, the probability of the kth label given an image (f) is

    P(ℳ_k| f) =P( f,ℳ_k)/P( f)=∫_ z P( f| z,ℳ_k) p( z|ℳ_k) d z P(ℳ_k)/∑_l=0^9 ∫_ z P( f| z,ℳ_l) p( z|ℳ_l) d z P(ℳ_l)


The above integrals are intractable, so are approximated.
This was done by a “Laplace approximation”[While many readers will be familiar with optimising the residual squared error (RSE), they may not understand how this differs from the evidence lower bound (ELBO). For a textbook explanation of Bayesian approaches, including the Laplace approximation, see <cit.>, <cit.> or <cit.>.] whereby the approximate distribution of z is given by

    q( z) =𝒩( z|ẑ, S^-1)

From this approximation, we can compute

    ∫_ z P( f,  z | ℳ) d z≃    P( f, ẑ | ℳ) ∫_ zexp(-12
    ( z-ẑ)^T  S ( z-ẑ) ) d z=    P( f, ẑ | ℳ)   | S/(2π)|^1/2


For each image (f), the mode (ẑ) of p( f, z | ℳ_k) was computed (see Section <ref>) by

    ẑ=_ z(
    J( f,  z, μ, Ŵ^a, Ŵ^v)
    + 12 z^T (λ_1 Â + λ_2 (Ŵ^a)^T  L^a Ŵ^a + λ_2 (Ŵ^v)^T  L^v Ŵ^v) z).

The Hessian of the objective function around this mode (Section <ref>) was used to approximate the uncertainty (S^-1).

Training was done with different sized subsets (300, 500, 1,000, 3,000, 5,000, 10,000, and all 60,000) of the MNIST training data, whereas testing was always done using the 10,000 test images.
In each of the training subsets, the first of the images were always used, which generally leads to slightly different sized training sets for each of the digits.
Example images, along with the fit from the models trained using the first 10,000 images, are shown in Fig. <ref>.
Model fitting was run for 20 iterations, using a Bernoulli likelihood with K=16, ν_0 = 16, λ = [0.95  0.05], ω^a = [0.002  0.2  0], ω^μ = N [10^-7 10^-5 0] and ω^v = [0.002  0.02  2  0.2  0.2].



When applied to medical images, machine learning can suffer from the curse of dimensionality.
The number of pixels or voxels in each image (M) is often much greater than the number of labelled images (N) available for training.
For MNIST, there are 60,000 training images, each containing 784 pixels, giving N/M ≃ 75.
In contrast, even after down-sampling to a lower resolution, a 3D MRI scan contains in the order of 20,000,000 voxels.
Achieving a similar N/M as for MNIST would require about 1.5 billion labelled images, which clearly is not feasible.
For this reason, this section focuses on classification methods trained using smaller subsets of the MNIST training data.
Accuracies are compared against those reported by <cit.> for their Deeply Supervised Nets, which is a deep learning approach that performs close to state-of-the-art (for 2015), particularly for smaller training sets.
Invariant scattering convolutional networks are also known to work well for smaller training sets, so some accuracies taken from <cit.> are also included in the comparison.
We are not aware of more recent papers that assess the accuracy of deep learning using smaller training sets.

Plots of error rate against training set size are shown in Fig. <ref>, along with the approximate error rates from <cit.> and <cit.>.
The plot shows the proposed method to be more accurate than deep learning for smaller training sets, but it is less accurate when using the full training set, as the error rate plateaus to a value of about 0.85% for training set sizes of around 5,000 onward.
Visual assessment of the fits to the misclassified digits (Fig. <ref>) suggests that relatively few of the failures can be attributed to registration errors.



These experiments with MNIST suggest that one avenue of further work could be to elaborate on the simple multivariate Gaussian model for the distribution of latent variables.
Although accuracies were relatively good for smaller training sets, the Gaussian assumptions meant that increasing the amount of training data beyond about 5,000 examples did not bring any additional accuracy.
One example of where the Gaussian distribution fails is when attempting to deal with sevens written either with or without a bar through them, which clearly requires some form of bimodal distribution to describe (see Fig. <ref>).
One approach to achieving a more flexible model of the latent variable probability density would to use a Gaussian Mixture Model (GMM) <cit.>.
This type of approach has been used for appearance models <cit.>.
In principle, a variational Bayesian GMM <cit.> could be used to automatically select the optimal model complexity, leading to an approach that self-tunes its complexity according to the amount and quality of of data available.
This type of model selection has previously been used for shape modelling <cit.>, as well as for other aspects of medical image computing.
With a more flexible model, it may be possible to achieve accuracies similar to those achieved by deep learning, but with fewer labelled training examples.
One of the aims of the Medical Informatics Platform of the HBP was to cluster patients into different sub-groups.
In addition to achieving greater accuracy (probably), incorporating a GMM over the latent variables could also lead to this clustering goal being achieved.






 §.§ Experiments with segmented MRI

Experiments were performed using 1,913 T1-weighted MR images from the following datasets.


  * The IXI dataset, which is available  under the Creative Commons CC BY-SA 3.0 license from <http://brain-development.org/ixi-dataset/>.
Information about scanner parameters and subject demographics are also available from the web site.  Scans were collected on three different scanners using a variety of MR sequences.
This work used only the 581 T1-weighted scans.

  * The OASIS Longitudinal dataset is described in <cit.>.  The dataset contains longitudinal T1-weighted MRI scans of elderly subjects, some of whom had dementia.  Only data from the first 82 subjects of this dataset were downloaded from <http://www.oasis-brains.org/>, and averages of the scans acquired at the first time point were used. 

  * The COBRE (Centre for Biomedical Research Excellence) dataset are available for download from <http://fcon_1000.projects.nitrc.org/indi/retro/cobre.html> under the Creative Commons CC BY-NC license. The dataset includes fMRI and T1-weighted scans of 72 patients with Schizophrenia and 74 healthy controls. Only the T1-weighted scans were used.  Information about scanner parameters and subject demographics is available from the web site.

  * The ABIDE I (Autism Brain Imaging Date Exchange) dataset was downloaded via <http://fcon_1000.projects.nitrc.org/indi/abide/abide_I.html> and is available under the Creative Commons CC BY-NC-SA license. There were scans from 1,102 subjects, where 531 were individuals on the Autism Spectrum. Subjects were drawn from a wide age range and were scanned at 17 different sites around the world.  All the T1-weighted scans were used, and these had a very wide range of image properties, resolutions and fields of view.  For example, many of the scans did not cover the cerebellum.


The images were segmented using the algorithm in SPM12, which uses the approach described in <cit.>, but with some additional modifications that are described in the appendices of <cit.>.
Binary maps of grey and white matter were approximately aligned into ICBM152 space using a rigid-body transform obtained from a weighted Procrustes analysis <cit.> of the deformations estimated by the segmentation algorithm.
These approximately aligned images have an isotropic resolution of 2 mm.




  §.§.§ 2D experiments with segmented MRI

It is generally easier to visualise how an algorithm is working when it is run in 2D, rather than 3D.
The examples here will be used to illustrate the behaviour of the algorithm under topological changes, when variability can not be modelled only via diffeomorphic deformations.

A single slice was extracted from the grey and white matter images of each of the 1,913 subjects, and the joint shape and appearance model was fit to the data using the settings for categorical image data.
This assumed that each voxel was a categorical variable indicating one of three tissue classes (grey and white matter, as well as background).  Each 2D image was encoded by 100 latent variables (i.e. K=100).
Eight iterations of the algorithm were used, with λ = [0.9  0.1], ω^a = [0.1  16  128], ω^μ = N [0.0001  0.01  0.1], ω^v = [0.001  0  32  0.25  0.5] and ν_0 = 100.



Some model fits are shown in Fig. <ref>, and the principal modes of variability are shown in Fig. <ref>.
As can be seen, these images are reasonably well modelled, although achieving a similar quality of fit for the full 3D data would probably require about 1,000 (100^3/2) variables.
Fitting a model of this size would require many more than the 1,900 subjects in this dataset.
Note that the topology of the images may differ, which (by definition[Topology is concerned with properties that are preserved following diffeomorphic deformations (see <https://en.wikipedia.org/wiki/Topology>).]) is not something that can be modelled by diffeomorphisms alone.
The inclusion of the appearance model allows these topology differences to be more accurately captured.








  §.§.§ Imputing missing data 

The ability to elegantly handle missing data is a useful requirement for mining hospital scans.
These often have limited fields of view, and may miss out parts of the brain that are present in other images.
The objective here is to demonstrate that a reasonable image factorisation can be learned, even when some images in the dataset may not have full organ coverage.

This experiment used the same slice through the data as above, and a rectangle covering 25% of the area of the images was placed randomly in each and every image of the training set (wrapping around at the edge of the field of view), and the intensities within these rectangles set to NaN (“not a number” in the IEEE 754 floating-point standard).
The algorithm was trained, using the same settings as described previously, on the these modified images.
Although imputed missing values may not be explicitly required, they do provide a useful illustration of how well the model works in less than ideal situations.
Fig. <ref> shows a selection of the images with regions set to NaN, and the same images with the missing values predicted by the algorithm.

The ability to handle missing data allows cross-validation to be used to determine the accuracy of a model, and how well it generalises.
In addition to the joint shape and appearance model, this work also allows simplified versions to be fitted that involve only shape (i.e., not using W^a, as in <cit.>) or in a form that varies only the appearance (i.e. not using W^v).
In addition, this work also includes a version where different sets of latent variables control the shape and appearance.
Here, there were 30 variables to control appearance, and 70 to control shape.
The aim was to compare the four models by assessing how well they are able to predict data that was unavailable to the model during fitting.
This gives us ground truth with which to compare the models' predictions, and is essentially a form of cross-validation procedure.
Accuracy was measured by the log-likelihood of the ground truth data, which was computed only for pixels that the models did not have access to during training.



The results of the cross-validation are shown in Fig. <ref>, and clearly show that the models that combine both shape and appearance have greater predictive validity than either the shape or appearance models alone.
Mean squared errors from the different models are also presented, and these exhibit the same general pattern.
Although the difference was small, the best results were from the model where each latent variable controls both shape and appearance, rather than when they are controlled separately.
Both of these combined models outperformed models of only shape or only appearance variability.
Note that changes to hyper-parameter settings, etc. may improve accuracies further.






  §.§.§ 3D experiments with segmented MRI

The aim of this section was to apply the method to a large set of 3D images, and use the resulting latent variables as features for pattern recognition.
For this, a version of the model was used whereby some latent variables controlled appearance, whereas others controlled shape.
The motivation for this was that it allows the different types of features to be differentially weighted when they are used to make predictions.

The algorithm was run on the full 3D dataset, using 70 variables to control shape and 30 to control appearance.  Eight iterations were used, with λ = [1  1], ω^a = [0.01  1  50], ω^μ = N [0.00001  0.01  0.1] and ω^v = [0.001  0  10  0.1  0.2].
The resulting model fits are shown in Fig. <ref>, as well as reconstructions using only the appearance part of the model or the shape model part.
This result can be compared against the 2D model fit shown in Fig. <ref>.
There are two reasons why the 3D fit explains a smaller proportion of the variability than for the 2D examples. The first is that the 3D model fit uses different variables to control shape and appearance, meaning that each can explain slightly less of the variability.
The second reason is simply that it is a 3D fit, so that there is a great deal more variability to explain.



The main objective of this work is to extract features from sets of medical images, which are effective for machine learning applications.
Here, a five-fold cross-validation is used to assess the effectiveness of these features.
Machine learning used a linear Gaussian process classification procedure, which is essentially equivalent to a Bayesian approach to logistic regression.
The implementation was based on the method for binary classification using expectation propagation described in <cit.>.
For the COBRE dataset, classification involved separating controls from patients with schizophrenia.
Similarly, the analysis of the ABIDE dataset involved identifying those subjects on the autism spectrum, with features orthogonalised with respect to the different sites.
Classification involved three hyper-parameters, which weighted the contributions from shape features, appearance features and a constant offset.

ROC curves are shown in Fig. <ref>.
For ABIDE, the accuracy was 57.6%.  While this is not especially high, it is close to the accuracy reported by others who have applied machine learning to the T1-weighted scans.
For example, <cit.> reported 60.1% classification accuracy with the same data using a support-vector machine with radial basis functions.
The accuracy achieved for the COBRE dataset was 74.7%, which is similar to the  69.7% accuracy reported by <cit.> using COBRE.  <cit.> achieved 71.4% accuracy for separating controls from subjects with schizophrenia, but using a different (larger) dataset of T1-weighted scans.
Both results were comparable to those obtained by <cit.>.









§ DISCUSSION

This work presents a very general generative framework that may have widespread use within the medical imaging community, particularly for those situations where conventional image registration approaches are more likely to fail.
As in many image analysis procedures, there is a certain amount of circularity arising from the dependencies among variables.
These circularities are not handled well by using a single pass through a conventional pipeline.
In contrast, by formulating the problem within a generative modelling framework, the algorithm is better able to deal with them.
Because of its generality, the model we presented should provide a good starting point for a number of avenues of further development, in addition to those previously mentioned that relate to going beyond a simple Gaussian distribution for the latent variables.

Most image analysis applications have a number of settings to be tuned, and the current approach is no exception.
Although this tuning is rarely discussed in papers, the settings can have quite a large impact on any results.
The effects of different forms of regularisation used for registration are illustrated in <cit.>.
The choice of settings plays a more important role when training using only a small number of images.
With more images, the contribution of these priors becomes less influential.
One aspect of the work that may need additional attention is the setting of λ_1 and λ_2.
The ideal settings would be λ_1=1 and λ_2=0, but in practice, greater regularisation is required in order to achieve good results.
A plausible explanation for this would be that assumptions of i.i.d. noise are not generally met, so a “virtual decimation factor”, which accounts for correlations among residuals, may need to be accounted for <cit.>.
The fact that the approach is not fully Bayesian (i.e., it only makes point estimates of many parameters and latent variables, rather than properly accounting for their uncertainty) is likely to be another reason why additional regularisation is needed.
Greater regularisation also seems to help avoid local optima, suggesting that further investigation into coarse-to-fine strategies may be warranted.

For computational reasons, the uncertainty with which the latent variables (z) are estimated is not treated in a properly Bayesian way.
In theory, more accurate results would be achieved by marginalising with respect to these latent variables.
This strategy is relatively straightforward for PCA with simple Gaussian noise assumptions <cit.>, but becomes more difficult for other types of data.
Monte-Carlo approaches <cit.> are currently too slow for routine use on large datasets, although a reasonable alternative might be to use variational Bayesian methods.
For example, <cit.> used a local variational approximation of the log partition function to do a type of generalised PCA of binary data.
<cit.> used a related approach for modelling the shape variability of binary images, although this work used a different local approximation.
Local variational approaches appear to be more difficult to extend to categorical data with more than two classes, and the image deformation part of the model would probably lead to many further difficulties for variational Bayesian methods.

There are a number of directions in which the current work could be extended.
One of the first areas we plan to investigate is to integrate the approach with current strategies for generating tissue probability maps <cit.>.
Rather than first segmenting images into different tissue classes and then fitting the proposed model to these, it would be more elegant to properly combine the approaches into a single unified generative model of the original images.
In principle, this strategy could lead to more robust segmentation because the ways in which tissue priors are allowed to vary are constrained to be more biologically plausible.

Another avenue is to allow some shape variability beyond what can be encoded by the first few eigenmodes.
For example, <cit.> combined the eigenmode representation with a model of additional shape variability, giving a framework that is conceptually related to that of <cit.>, as this allows a covariance matrix over velocity fields to be defined and optimised.
Approaches building on <cit.>, who proposed to learn the residual covariance by assuming finite support (and thus sparsity) of the underlying covariance function, might also be potentially interesting.

The framework would also generalise further for handling paired or multi-view data, which could add a degree of supervision to the method.
There have been a number of publications on generating age- or gender-specific templates, or on geodesic regression approaches <cit.> for modelling trajectories of ageing.
Concepts from joint matrix factorisation approaches, such as canonical correlation analysis <cit.>, could be integrated into the current work, and these could be used to allow the model fitting to be informed by age, gender, disease status etc.




  §.§.§ Acknowledgements

This project has received funding from the European Union’s Horizon 2020 Research and Innovation Programme under Grant Agreement No. 720270 (HBP SGA1).
The Wellcome Centre for Human Neuroimaging is supported by core funding from the Wellcome Trust [grant number 203147/Z/16/Z].  Funding for the OASIS dataset came from the following grants: P50 AG05681, P01 AG03991, R01 AG021910, P20 MH071616, U24 RR021382.
Funding for the ABIDE I dataset came from a variety of sources, which include NIMH (K23MH087770 and R03MH096321), the Leon Levy Foundation, Joseph P. Healy and the Stavros Niarchos Foundation to the Child Mind Institute.







§ ALGORITHM 

A highly simplified version of the factorisation algorithm is shown in Algorithm <ref>.
It is an expectation maximisation approach, which involves alternating between computing the shape and appearance basis functions (plus a few other variables – M-step), and computing the expectations of the latent variables (E-step).
For better convergence of the basis function updates of the M-step, an orthogonalisation step is included in the algorithm.

The M-step relies on Gauss-Newton updates of three elements: the mean template (μ), shape subspace (W^a) and appearance subspace (W^v).
These updates have the general form of w w - ( H +  L)^-1( g +  L w), where L is a very sparse Toeplitz or circulant matrix encoding spatial regularisation, and H encodes a field of small matrices that are easy to invert.
The full-multigrid method, described in <cit.>, is particularly well suited to solving this type of problem.

The E-step involves updating the distributions of the latent variables (Z) and Gaussian prior (A).
To break the initial symmetry, the latent variables are all initialised randomly, while ensuring that ẐẐ^T is orthonormal.
This is easier than random initialisation of Ŵ^a and Ŵ^v, which are instead initialised to zero. 

In most of this appendix, matrices are written in bold upper-case (e.g., W^a, Z, etc).
In the computations, images are treated as vectors. These are written as lower-case bold, which includes the notation for individual columns of various matrices (e.g., w^a_k denotes the kth column of W^a, z_n denotes the nth column of Z, etc).
Scalars are written in italic, with dimensions in upper-case.
Estimates or expectations of parameters are written with a circumflex (e.g., Ẑ).
Collections of vectors may be conceptualised as matrices, so are written in bold-upper-case (e.g., G^a, where individual vectors are g^a_k).
Collections of matrices are written in “mathcal” font (e.g., ℋ^a, where individual matrices are H^a_kk).
There is no systematic rule about the use of Greek letters, although sometimes warping an image may be written as a(ψ) and at other times, warping may be conceptualised as matrix multiplication and written as Ψ a.
The matrix transpose operation is denoted by the “^T” superscript (as in Ψ^T).
Creating a diagonal matrix from a vector (as in (exp q)), as well as treating the diagonal elements of a matrix as a vector (as in ( Q)) are both denoted by “”. 
The trace of a matrix (sum of diagonal elements) is denoted by “”.
Sometimes, gradients of an image are required.  In 3D, the three components of the spatial gradient of a would be denoted by ∇_1  a, ∇_2  a and ∇_3  a.

Comments in Algorithm <ref> saying “Dist” indicate which steps should be modified for running within a distributed privacy-preserving framework.
The idea here is that the main procedure would be run on the “master” computer, whereas various functions would be run on the “worker” machines on which the data reside.
These workers would only pass aggregate data back to the master, whereas the latent variables, which explicitly encode information about individuals, would remain on the workers.
As the algorithm is described here, the images (F) and expectations of the latent variables Ẑ are passed back and forth between the master and workers, but this need not be the case.
If these data and variables were all to reside on the worker machines, the master machine would still be able to run using only the aggregate data.

For simplicity, Algorithm <ref> does not include functions for computing variances (σ^2 used by the Gaussian noise model), etc., and these variables are not shown to be passed to the various functions that use them.
However, it should be easy to see how these changes would be incorporated in practice.

Also, the illustration does not show any steps requiring the objective function, which include various backtracking line-searches to ensure that parameter updates cause the objective function to improve each time.
In practice, the algorithm is run for a fixed number of iterations, although the log-likelihood could be used to determine when to stop.







 §.§ Updating the mean (μ̂)

From Eqn. <ref>, we see that a point estimate of the mean (μ) may be computed by

    μ̂=_μ(
        12μ^T  L^μμ
        + ∑_n=1^N J( f_n, ẑ_n, μ, Ŵ^a, Ŵ^v)
    ).


In practice, this log probability is not fully maximised with respect to μ at each iteration.
Instead, μ̂ is updated by a single Gauss-Newton iteration.
This requires gradients and Hessians computed as shown in Algorithm <ref>, which simply involves summing over those computed for the individual images.
A small amount of regularisation is used for the estimate of the mean.
One of the reasons for using this is that it alleviates some of the problems with the gradients of Jacobian-weighted averaging that was pointed out in Appendix B of <cit.>, thus leading to better convergence.
This is especially important in situations where it can help to smooth over some of the effects of missing data.






 §.§ Likelihood derivatives

The algorithm can be run using a number of different appearance models, and the gradients and Hessians involved in the Gauss-Newton updates depend upon the one used.




  §.§.§ Gaussian model

Algorithm <ref> shows derivatives for the Gaussian noise model.
For a single voxel, this is based on

    J_L_2   =12ln(2π) + 12lnσ^2 + 12σ^2||f - a'||_2^2 
    d J_L_2/d a'   =1σ^2 (a' - f)  and d^2 J_L_2/d a'^2 = 1σ^2

For voxels where data is missing, both J_L_2 and d J_L_2d a' are assumed to be zero.
Using matrix notation, the objective function for an image is therefore

    J' =12 σ^2 (Ψ a -  f)^T(Ψ a -  f) + M2(ln(σ^2) + ln(2π)).

The gradients and Hessians, with respect to variations in a, are

    g' =   Ψ^T (1σ^2 (Ψ a -  f) )
     H' =   1σ^2Ψ^T Ψ

In practice, the Hessian (H') is approximated by a diagonal matrix

    H' ≃   ( Ψ^T  11σ^2)

where 1 is a vector of ones.
This approximation works in the optimisation because all rows of Ψ sum to 1, so for any vector d of the right dimension, the rows of Ψ^T ( d) Ψ sum to Ψ^T  d.
Because (for trilinear interpolation) all elements of Ψ are greater than or equal to zero, so if all elements of d are non-negative, then all eigenvalues of ( Ψ^T  d) - Ψ^T ( d) Ψ are greater than or equal to zero.
Non-negative eigenvalues ensure that the approximation to the Hessian is positive semi-definite, and thus the optimisation aims for a minimum, rather than a maximum.






  §.§.§ Binary model

Slight modifications are made to the algorithm in order to use the Bernoulli noise model with the sigmoidal squashing function.
For each voxel, the negative log-likelihood is

    J_Bern   = -(f a' + ln s(-a')) , where  s(a') = 1/1+exp(-a').

Modifications are made to the gradient and Hessian of Algorithm <ref>, based on the derivatives

    d J_Bern/d a'   = s(a') - f  and d^2 J_Bern/d a'^2= s(a')(1-s(a')).

Using matrix notation (where s≡ s( a)), the gradients and Hessians are

    g' =   Ψ^T (Ψ s -  f)
     H' =   Ψ^T ( s)(1- s) Ψ≃( Ψ^T ( s) (1- s) )





  §.§.§ Categorical model

For the categorical model with a softmax squashing function, the negative log-likelihood of a single voxel is

    J_cat   = -∑_k=1^K a'_k f_k + log(∑_k=1^K exp a'_k)


This would use the gradients and Hessians

    d J_cat/d a'_k   = s_k( a') - f_k , where  s( a') = exp a')/∑_k=1^K exp a'_k
    d^2 J_cat/d a'_k a'_l   = s_k( a') (δ_jk - s_j( a'))

where δ_jk is the Kronecker delta function.
Computation of the gradients and the approximation of the Hessian follow similar lines to those for the binary and Gaussian models. 




 §.§ Geodesic shooting

The algorithm requires a means of computing diffeomorphic deformations from the initial velocities via a Geodesic shooting procedure.  Algorithm <ref> shows how this is achieved.
In the presented algorithm, Dψ denotes the Jacobian tensor field of ψ, and (Dψ)^T u indicates a pointwise multiplication with the transpose of the Jacobian.
|Dψ| denotes the field of Jacobian determinants.
L v in the continuous framework is equivalent to the matrix multiplication L^v  v in the discrete framework. 
The operation L^g u denotes applying the inverse of L to u, such that L L^g u = u.
In practice, this is a deconvolution, which is computed using fast Fourier transforms.
Much has already been written about the geodesic shooting procedure, so the reader is referred to <cit.> for further information.






 §.§ Updating appearance basis functions (Ŵ^a)

Appearance basis functions are optimised by

    Ŵ^a =_ W^a(
        12((λ_1 N I + λ_2 ẐẐ^T) ( W^a)^T  L^a  W^a )
        + ∑_n=1^N J( f_n, ẑ_n, μ̂,  W^a, Ŵ^v)
    ).


The first step involves computing the gradients and Hessians, which can be done in a distributed way and is shown in Algorithm <ref>.
Note that Algorithm <ref> only shows the computation of gradients and Hessians for the Gaussian noise model, and that slight modifications are required when using other forms of appearance model.
Gradients and Hessians for updating these basis functions (W^a) are similar to those for the mean updates, except for weighting based on the current estimates of the latent variables.
Note that for this approach to work effectively, the rows of Ẑ should be orthogonal to each other, which is explained further in Section <ref>.
Note that only a single Gauss-Newton iteration is performed, so the objective function in Eqn. <ref> is not fully optimised, but merely improved over its previous value. 







 §.§ Updating shape basis functions (Ŵ^v)

Shape basis functions are optimised by

    Ŵ^v =_ W^v(
        12((λ_1 N  I + λ_2 ẐẐ^T) ( W^v)^T  L^v  W^v )
        + ∑_n=1^N J( f_n, ẑ_n, μ̂, Ŵ^a,  W^v)
    ).


A single Gauss-Newton iteration is used to update the basis functions of the shape model (W^v), which is done in such a way that changes to W^v improve the objective function with respect to its previous value, rather than fully optimise Eqn. <ref>.
While most Gauss-Newton iterations improve the fit, a backtracking line search is included to ensure that they do not overshoot.
For simplicity, this aspect of the algorithm is not shown in Algorithm <ref>.
As for updating W^a, this requires the rows of Ẑ to be orthogonal to each other.

As for updating the appearance basis functions, computing the gradients and Hessians for the shape bases can be also done in a distributed way.
The strategy for computing gradients and Hessians is shown in Algorithm <ref>.






 §.§ Updating latent variables (ẑ_n) 

Within a distributed multi-centre privacy-preserving framework, the strategy would be to have the basis functions shared across sites, whereas the image data and latent variables would remain hidden within each of the sites.
All computations that use the image data and latent variables would be done at each site, with only aggregate data shared across them.
The modes of the latent variables are updated via a Gauss-Newton scheme, similar to that used by <cit.>.


    ẑ_n =_ z_n(
    J( f_n,  z_n, μ, Ŵ^a, Ŵ^v)
    + 12 z_n^T (λ_1 Â + λ_2 (Ŵ^a)^T  L^a Ŵ^a + λ_2 (Ŵ^v)^T  L^v Ŵ^v) z_n
    )




The inverse of the (approximate) Hessians allows a Gaussian approximation of the uncertainty with which the latent variables are updated to be computed (“Laplace approximation”).
This is the S matrix (returned by Algorithm <ref>), which is combined with ẐẐ^T (returned as C^z) and used to re-compute Â.




 §.§ Expectation of the precision matrix (Â) 

This work uses a variational Bayesian approach for approximating the distribution of A, which is a method described in more detail by textbooks, such as <cit.> or <cit.>.
Briefly, it involves taking the joint probability of Eqn. <ref>, discarding terms that do not involve A, and substituting the expectations of the other parameters into the expression.
This leads to the following approximating distribution, which can be recognised as Wishart.

    ln q( A) =   12(N + ν_0 - K - 1) ln| A| - 12(([ Z Z^T] + Λ_0^-1)  A) + const=   ln𝒲( A | Λ, ν)

where Λ = ([ Z Z^T] + Λ_0^-1)^-1  and ν = ν_0 + N.
In practice, [ Z Z^T] is approximated by C^z +  S, described previously.

Other steps in the algorithm use the expectation of A, which (see Appendix B of <cit.>) is

    Â=[ A] =νΛ.





 §.§ Orthogonalisation 

The strategy for updating Ŵ^a and Ŵ^v involves some approximations, which are needed in order to save memory and computation. This approximation is related to the Jacobi iterative method for determining the solutions to linear equations, which is only guaranteed to converge for diagonally dominant matrices.
Rather than work with the Hessian for the entire W matrix together, only the Hessians for each column of W are computed by Algorithms <ref> and <ref>.
This corresponds with a block diagonal Hessian matrix for the entire W, which has the form

    H=[ H_11    0       0 0 H_22        0⋮    ⋮    ⋱  ⋮ 0    0      H_KK ].


More stable convergence can be achieved by transforming the basis functions and latent variables in order to minimise the amount of signal that would be in the off-diagonal blocks, thus increasing the diagonal dominance of the system of equations.
In situations where diagonal dominance is violated, convergence can still be achieved by decreasing the update step size.
This is analogous to using a weighted Jacobi iteration, where in practice the weights are found using a backtracking line-search.

Signal in the off-diagonal blocks is reduced by orthogonalising the rows of Ẑ.
This is achieved by finding a transformation, T, such that TẐ ( TẐ)^T and  (Ŵ^v T^-1)^T  L^v Ŵ^v T^-1 + (Ŵ^a T^-1)^T  L^a Ŵ^a T^-1 are both diagonal matrices.
Transformation T is derived from an eigendecomposition of the sufficient statistics, whereby the symmetric positive definite matrices are decomposed into diagonal (D^z and D^w) and orthonormal (V^z and V^w) matrices, such that

    V^z  D^z ( V^z)^T =    C^z,
     V^w  D^w ( V^w)^T =    C,

where C^z = ẐẐ^T and C = (Ŵ^v)^T  L^v Ŵ^v + (Ŵ^a)^T  L^a Ŵ^a.

A further singular value decomposition is then used, giving

    U D V^T = ( D^w)^1/2 ( V^w)^T  V^z ( D^z)^1/2.


The combination of various matrices is used to give an initial estimate of the transform

    T= D V^T ( D^z)^-1/2 ( V^z)^T .


The above T matrix could be used to render the matrices orthogonal, but their relative scalings would not be optimal. 
The remainder of the orthogonalisation procedure involves an iterative strategy similar to expectation maximisation, where the aim is to estimate some diagonal scaling matrix Q with which to multiply T.
This matrix is parameterised by a set of parameters q, such that

    Q=(exp q).


The first step of the iterative scheme involves re-computing Â, as described in Section <ref>, but incorporating the current estimates of Q T.

    Â=νΛ =  (N + ν_0)( Q T ( C^z +  S)( Q T)^T + Λ_0^-1)^-1.


The next step in the iterative scheme is to re-estimate q, such that

    q̂=_ q    ( ( (exp (- q)) ( T^-1)^T C T^-1(exp (- q)) ) 
                 +    ((exp q) T C^z  T^T(exp q) Â) ).

This can be achieved via a Gauss-Newton update, which uses first and second derivatives with respect to q.
The overall strategy is illustrated in Algorithm <ref>, which empirically is found to converge well.




