Bayesian shape modelling of cross-sectional geological data
Thomai Tsiftsi Ian H. Jermyn Jochen Einbeck








Bayesian shape modelling



Tsiftsi et al.    
                                    
                                    



University of Durham, Department of Mathematical Sciences, Statistics and Probability Group, UK


thomai.tsiftsi@durham.ac.uk




Shape information is of great importance in many applications. For example, the oil-bearing capacity of sand bodies, the subterranean remnants of ancient rivers, is related to their cross-sectional shapes. The analysis of these shapes is therefore of some interest, but current classifications are simplistic and ad hoc. In this paper, we describe the first steps towards a coherent statistical analysis of these shapes by deriving the integrated likelihood for data shapes given class parameters. The result is of interest beyond this particular application. 





Bayesian shape modelling of cross-sectional geological data
    Thomai Tsiftsi1, Ian H. Jermyn1, Jochen Einbeck1
    December 30, 2023
===========================================================









§ INTRODUCTION

Sand bodies, the sedimentary, subterranean remnants of ancient rivers, are important to both geology and the petroleum industry. In particular, their cross-sectional shapes help determine their oil-bearing capacity. Current classification schemes for sand body shapes are qualitative, simple, and ad hoc, and so there is a need for a quantitative analysis with the help of statistical models. There are several problems of interest: estimation of shape class parameters given labelled data shapes (a `data shape' is an ordered set of points in ℝ^2); classification of new data shapes; and unsupervised classification. Parameter estimation is described by the probability ℙ(w|y, c), where w denotes the shape class parameters and y the dataset, which consists of several data shapes, together with their class labels c. By Bayes' theorem, this is given by: 
    
    ℙ(w|y, c)∝ℙ(y|w, c)ℙ(w)  .


In this, as in all of the above problems, the major task is to calculate the likelihood ℙ(y|w, c). This is the problem addressed in this paper. The problem is not unique to the sand body application: it occurs in many applications of shape modelling, and is thus of broad interest. 




§ THE LIKELIHOOD

The calculation of the likelihood is complicated due to the presence of many nuisance parameters that must be integrated over. The partitioned likelihood is: 
    
    ℙ(y|w, c) = 
    ∑_b ∈ B∫ ds  dg  dβ dσ ℙ(y|σ, b, s, g, β)
            ℙ(σ) 
            ℙ(b) 
            ℙ(s|β)
            ℙ(g)  
            ℙ(β|w, c)  ,


where we have made a number of simplifying independence assumptions. Here  ℙ(y|σ, b, s, g, β) = exp(-1/2σ^2∑_i=1^N |y(i)-g∘β (s(b^-1(i)))|^2 ) models errors in shape point collection as Gaussian white noise.

In the above expression, β is the underlying sand body shape modulo similarity transformations, which comes from a class C with parameters w, while g ∈ G ≡ SO(2)⋉𝐑^2 ×𝐑_+ is a similarity transformation generating the full sand body shape gβ. Data formation is modelled as a sampling s of N points around the sand body shape, and a bijection B ∋ b: [1,...n ]→ [1,...n] relating each point of the sand body shape to a unique point of the data shape, giving gβ(s(b^-1)), plus the above Gaussian noise with variance σ.

In previous work, e.g. Dryden and Mardia (1998) or Srivastava and Jermyn (2009), an algorithmic approach was taken to the integrals over the group G, using the Procrustes algorithm to compute a zeroth order Laplace approximation. Here we carry out the group integrations, and the integration over σ, analytically, resulting in a closed form expression. This is the main contribution of this paper.

First, we have to choose the priors for g and σ. Jeffrey's joint prior for g and σ was calculated to be a Var(β(s))/σ^5, but this leads to a non-normalizable likelihood. Instead, a regularized version was employed. A Gaussian prior was used for translations; a uniform prior for rotation angle; and a Rayleigh prior for scaling. These choices break translation invariance by effectively limiting the size of the two-dimensional domain in ℝ^2 in which the shape points lie, and break scale invariance by effectively limiting the range of scales considered. With these priors, the result of the integrations over translations, rotations, and scalings is: 
    
    ℙ(y|b, s, β) = 
    1/Z1/σ^2nexp{-1/2σ^2(
                    ñṼãr̃(̃ỹ)̃
                    -ñ^2C̃õṽ^̃2̃(β                 (s(b^-1)),y)/ñṼãr̃                 (β (s(b^-1)) + 1/B^2)
            } ,


where:


  * 

ñ=n+1/D^2


  * 

Ṽãr̃(y)=1/ñ[∑_i|y_i|^2-1/ñ∑_i∑_jy_iy_j]


  * 

C̃õṽ(β (s(b^-1)),y)=1/ñ[ ∑_iβ (s(b^-1(i)))y_i-1/ñ∑_i∑_jβ (s(b^-1(i))y_j]



which are regularized versions of the number of points and the variance. B, D, α, c are appropriate regulators and Z the normalization constant. 

A Γ prior was used for σ. After integration, this leads to:
    
    ℙ(y|b, s, β) = 
            Γ(n+α)/Z×[  
                ñṼãr̃(̃ỹ)̃-ñ^2             |C̃õṽ(β (s(b^-1)),y)|^2/ñṼãr̃(β) 
                (s(b^-1)) + 1/B^2 + 2c   
            ]^n-α .


This expression is the main result of the paper. It applies to any shape modelling application in which white Gaussian noise is added to a discrete set of shape points. 

To finish the calculation of the likelihood, we have to perform the s and β integrations, and the b summation, in equation (<ref>). This we do using Monte Carlo techniques. We use a uniform distribution on [0, 1]^N for samplings s, while β, which consists of positive quantities such as aspect ratios and lengths, is modelled with Γ distributions, whose parameters over all classes constitute w. In previous work, when the group integrations were approximated using the Laplace approximation and the Procrustes algorithm, the sum over bijections could be approximated, again in a zeroth order Laplace estimation, using the Hungarian algorithm, since only a linear assignment problem was involved. In the integrated likelihood derived here, the terms involving C̃õṽ(̃β̃ ̃(̃s̃(̃b̃^̃-̃1̃)̃)̃,̃ỹ)̃ complicate the situation, and turn the linear assignment problem into quadratic assignment, which is NP-hard. Instead of using the Laplace approximation, we approximate the full summation using Monte Carlo, with a uniform prior on b. 



§ PARAMETER ESTIMATION

The above result can be used to estimate the class parameters w given data shapes from each class. Figure <ref> shows an example of a likelihood surface, in this case computed for the simplest case of a rectangle modulo similarities. This is a one-dimensional shape space that can be parametrized by the aspect ratio. There is thus one Gamma distribution involved, and w is two-dimensional. The surface is rough, since only a coarse grid in parameter space was used to reduce simulation time. 

Maximum likelihood estimation was carried out using a built-in Matlab optimization function. The algorithm correctly navigated towards the maximum of the likelihood surface, but convergence was slow. We hope to improve this in the future by using a different maximization algorithm. 




§ CONCLUSION

The main contribution of this paper is the analytical evaluation of the integrals over the similarity group and the noise variance in a model for shape data. The application considered here is to the classification of the cross-sectional shapes of sand bodies, but the same techniques apply to any shape model involving Gaussian noise. This work is still in progress but shows promising signs of improving the current classification and estimation methods employed by geologists. There are technical obstacles in the form of numerical integrations which we hope to overcome in the near future.



  Dryden, I.L. and Mardia, K (1998).
    Statistical shape analysis.
    J. Wiley.

  Srivastava, A. and Jermyn, I.H. (2009)
    Looking for shapes in two-dimensional, cluttered point clouds.
    IEEE Trans. Patt. Anal. Mach. Intell., 31(9),
    1616 – 1629.



Notes 

In the manuscript version published in the conference proceedings volume, there was a small error in expression (4). This is now corrected.


