It is assumed that an initial correspondence field D, that matches the source shapeSsonto the target shapeStexists, wherexiis a point on the source shape and thereforexi+di is the corresponding point on the target surface. This can be created for each vertex in the source mesh by finding the closest point on the target shape. In section C.4.1, it is explained in details how an initial cor-respondence field is found. The problem of finding the optimal corcor-respondence vector field is then cast into a Bayesian framework of Markov Random Field restoration.
C.3.1 Bayesian Paradigm
The Bayesian paradigm proposed by Besag and others consists of four successive stages [24]. Adapted to the current application they become:
C.3 Markov Random Field Regularisation of Correspondences 127 1: Construction of a prior probability distribution p(d) for the
correspon-dence fieldD.
2: Formulation of an observation modelp(y|d) that describes the distribution of the observed shapes Y given any particular realization of the prior distribution. Y is the combined source and target shape.
3: Combination of the prior and the observation model into the posterior distribution by Bayes theorem
p(d|y) = p(y|d)p(d)
p(y) . (C.1)
4: Drawing inference based on the posterior distribution.
The goal is to find the correspondence field that maximises the posterior prob-ability
dˆ = arg max
d p(d|y). (C.2)
In the following, it is explained how the above-mentioned probability distribu-tions are defined using the theory of Markov Random Fields.
C.3.2 Graphs and Neighbourhoods
To describe a probability distribution on a spatial arrangement of points some useful definitions from graph theory are necessary.
Given a graph of n connected sites S = {si}ni=1. A neighbourhood system N ={Ns, s∈ S} is any collection of subsets ofS for which i) s /∈Ns, and ii) r ∈Ns⇔s ∈Nr, then Ns are the neighbours of s. A clique C is a subset of sitesS for which every pair of sites are neighbours. The surface mesh is treated as an undirected graph and i ∼ j is used to denote that sitei and sitej are neighbours.
C.3.3 Gibbs Distributions and Markov Random Fields
Given a neighbourhood systemN on the set of sitesS, the probability distrib-ution of any family of random variables indexed byS, i.e. D={Ds|s∈S}can now be considered. In this paperD is the correspondence field andDsis a 3D vector placed at the vertex with indexs.
The Markov property of MRFs is expressed as [107]
p(ds|dr, r6=s) =p(ds|dr, r∈Ns)∀s∈S, (C.3) which means that the probability distribution of each site only depends on the state of neighbouring sites. This does not exclude long-range correlations in the probability distribution over the entire graph though. To be able to specify the conditional probabilities in practice the equivalence between the Gibbs distributions and Markov Random Fields is used.
Given a neighbourhood system N ={Ns} let all cliques be denoted by C. For allC∈ C it is assumed that a family of potential functionsVCexist. An energy function of any given configuration ofD can now be defined
U(d) =X
C∈C
Vc. (C.4)
This leads to the definition of the Gibbs distribution induced by the energy functionU(d)
p(d) = 1
Zexp(−U(d)/T), (C.5)
where Z is the partition function and T is a parameter referred to as tem-perature. The Gibbs distribution maximises entropy (uncertainty) among all distributions with the same expected energy. The temperature controls the
“peaking” of the density function. The normalising constant, Z, may be im-possible to obtain due to the curse of dimensionality but often only ratios of probabilities are needed and the constant cancels out. The relation between the Gibbs distribution and Markov Random Fields are specified in the following theorem
Theorem C.1 (Hammersley-Clifford). LetN be a neighbourhood system. Then D is a Markov random field with respect to N if and only if p(d) is a Gibbs distribution with respect toN.
A proof can be found in [106]. Thus the task is to specify potentials that induce the Gibbs distribution in order to encompass MRF properties ofD.
Since the goal is to model the correspondence between Ss and St the corre-spondence vectors are bound to the surfaces, in effect only posing constraints on the length of the vectors at the individual sites. In practice, the constraint may be enforced by projection of the end of the correspondence vectors onto the closest point of the target surface in every site update of the relaxation. This is explained in detail in section C.4.2.
C.3 Markov Random Field Regularisation of Correspondences 129
C.3.4 Prior Distributions
Similar to pixel priors [24] energy functions based on differences between neigh-bouring sites are constructed. Extending to the multivariate case the general expression of the energy governing the site-priors becomes
Usite(d) =X
i∼j
||di−dj||pp, (C.6)
where || · ||p is the p-norm, 1 ≤p ≤ 2, and di represents the correspondence vector at theith site.
Withp= 2 the energy function induces a Gaussian prior on the correspondence field as used by for example Besag [24] in image restoration. Neglecting regions with strong surface dynamics the local optimisation becomes convex and the maximum likelihood (ML) estimate of the correspondence vector at theith site is taken as the mean of the neighbouring correspondence vectors. Since the edges of the meshes have not unit lengths, it is necessary to use a weighted average
dˆi= X
j∈Ni
wijdj/ X
j∈Ni
wij. (C.7)
Herewij are Gaussian weights derived from a fixed kernel size and calculated on the basis of the Euclidean distance between vertexiandjin the source mesh. A pure prior MRF restoration of correspondences hence solves the aperture prob-lem and the three-dimensional interpolation probprob-lem simultaneously by finding the simplest correspondence field.
Using this prior the maximum a posteriori (MAP) estimate of the MRF opti-mised purely by Iterative Conditional Modes is similar to the steady state of the algorithm for geometry constrained diffusion (GCD) proposed by Andresen and Nielsen [8]. GCD is a numerical scheme for solving a space and time discretised version of the heat equation on the correspondence field with certain boundary conditions. The GCD implementation works on volume-voxel diffusion, while the MRF method described in this paper works with 2D surfaces embedded in a 3D space. Moreover, the GCD approach in essence applies gradient descend optimisation and is directly dependable on a good initialisation of the corre-spondence field. The MRF formulation provides a natural framework allowing for more advanced optimisation.
Abandoning homogeneity and isotropy of the MRF non-global kernels may be introduced. Thus, adaptive Gaussian smoothing may be applied, e.g. by setting the standard deviation of the kernel to the square-root of the edge length of the closest neighbour of site i on the graph. Moreover, using the p= 1 norm
induces a median prior, with the ML estimate being the median of the corre-spondence vectors at the weighted neighbouring sites. This property makes the MRF attractive for correspondence fields with discontinuities, thus avoiding the smearing of edges attained by the Gaussian prior.
C.3.5 Observation Models
The observation model p(y|d) describes the conditional distribution of the ob-served shapes Y. A mapping that makes correspondences between regions of similar surface properties may be favoured by specifying an observation model.
The similarity measures may include derived features of the observed surfaces such as the curvature, orientation of the surface normals, or even surface texture.
The energy functionUmodel must be specified and depends on the application but typically involves a linear combination of multiple derived measures. In the following a selection of surface features is proposed.
C.3.5.1 Local Surface Orientation
The simple dot product between the surface normals may form the basis in specifying a governing energy function that favours correspondence between regions of similar orientation by
Unorm(y|d) =X
i
||nTs,int,i−1||q, (C.8) wherens,iis the surface normal at vertexiwith positionxion the source surface Ss, andnt,i the normal of the target surfaceSt at the coordinatexi+di. The parameterq >0 controls the sensitivity of the energy function. The normal at xi+di is found by weighting the normals from the triangle in which xi+di
is located. This is done by using the points barycentric coordinates as done in for example Phong shading [258, 99]. This energy function is expected to be useful when the source and target shape are roughly pre-aligned since it is not rotation invariant. Results for this energy function can be found in [199].
C.3.5.2 Surface Curvature
The use of second order geometric invariants such as curvature allows matching of homologous surface points despite a slight rotation of the surface. There are
C.3 Markov Random Field Regularisation of Correspondences 131 at least three types of surface curvatures that can be computed at each vertex:
mean, Gaussian and total curvatures [82]. A general setting for the curvature energy is :
Ucurv(y|d) =X
i
||cs,i−ct,i||pp, (C.9) wherecs,i is a vector containing the differential properties of the source surface at locationxi and ct,i contains the differential properties of the target surface at the coordinate xi+di. The differential properties at xi+di is found by interpolating as in section C.3.5.1.
C.3.5.3 Surface Texture
When surface texture is available multivariate locally derived features on the source and target shapes, ts and tt, may be applied in an energy function operating in the p-norm
Utext(y|d) =X
i
||ts,i−tt,i||pp, (C.10) wherets,i is a vector containing the texture of the source surface at locationxi
andtt,i contains the texture of the target surface at the coordinatexi+di.
C.3.6 Posterior Probability Distribution
Normalisation of the energy terms from the different prior and observation mod-els is typically chosen such that they operate on the same domain. However, the data analyst may choose to favour some terms over others, e.g. by relaxing the smoothness conditions in favour of correspondences between regions with similar surface characteristics. The total energy is:
Utotal= (1−α)Umodel+αUsite, (C.11) in which α ∈ [0 : 1] weights the influence of the model terms and Umodel is a linear combination of measures as explained in section C.3.5. By combining equation (C.1), (C.5) and (C.11) the posteriori conditional probability distrib-ution is given by
p(d|y)∝exp(−Utotal/T). (C.12) The method used to find the MAP estimate
dˆ = arg max
d p(d|y), (C.13)
is dependent on the complexity of the energy function. In section C.4.2 it is explained how The MAP estimate of the correspondence field can be found using stochastic optimisation.