Thesis Title: Bayesian Computational Methods for Spatial Analysis of Images.
Candidate: Matthew T. Moores
Discipline: Statistical Science
Prof. Ian Turner
Prof. Tony Pettitt
Dr. Fiona Harden
Dr. Chris Drovandi
This thesis is motivated by an important applied problem in image-guided radiation therapy. The aim is to assist in interpreting cone-beam computed tomography (CT) of radiotherapy patients by labelling the image pixels according to tissue type. These medical images have poor contrast-to-noise ratio (CNR), thus external sources of information are needed for accurate segmentation. Such sources include the individualised treatment plan, which is based on a diagnostic-quality CT scan that has been manually labelled by a clinician. We also use published studies of physiological variability to derive an estimate of spatial uncertainty.
To address this problem, we adopt a hidden Potts model of the image lattice. We introduce a method for deriving and representing the spatial prior for each patient as an external field in the hidden Potts model. Tissue density estimates are derived from the planning CT and adjusted to account for differences in image modality, forming priors for the noise parameters. These priors can be updated sequentially as more images of the patient are acquired.
Scalable computational algorithms are required for Bayesian inference on 3D volumetric images of this size. We evaluate the existing methods for intractable likelihood problems, including path sampling, pseudolikelihood and the approximate exchange algorithm. We introduce a precomputation step that involves fitting a binding function between the parameters and sufficient statistic of the Potts model. Using this precomputation, we achieve two orders of magnitude improvement in the scalability of approximate Bayesian computation with sequential Monte Carlo (ABC-SMC).
Matthew completed a Bachelor of Information Technology at QUT in 1996 and a Master of Mathematical Science in 2008. Prior to commencing his PhD, he was involved with the Visible Cell project at the Institute for Molecular Bioscience, UQ. He has also worked in R&D for various international companies, including DEC, Compaq and Oracle. He was a member of the development teams for the FaceWorks facial animation software, SpeechBot multimedia search engine, AgileTV speech-enabled programme guide, and the Healthcare Transaction Base. Matthew’s research interests include computational statistics, Bayesian inference, and image analysis
I’ve been invited to give two seminars this week at UNSW. Full details available here:
1:30pm, Thursday July 3: Intro to approximate Bayesian computation (ABC)
4:10pm, Friday July 4: Pre-processing for ABC-SMC with undirected graphical models
ABC is a useful method for Bayesian inference with intractable likelihoods, but existing algorithms require a large number of simulations of pseudo-data from the generative model. The computational cost of these simulations can be prohibitive for high-dimensional data. In this talk I will show that the scalability of ABC-SMC can be improved by performing a pre-computation step before model fitting. The output of this pre-computation can be reused across multiple datasets. I illustrate the approach using undirected graphical models, including the hidden Potts model for image analysis of large datasets.
This is the first in a series of posts describing the functions and algorithms that I have implemented in the
R package bayesImageS.
Gibbs sampling was originally designed by Geman & Geman (1984) for drawing updates from the Gibbs distribution, hence the name. However, single-site Gibbs sampling exhibits poor mixing due to the posterior correlation between the pixel labels. Thus it is very slow to converge when the correlation (controlled by the inverse temperature β) is high.
The algorithm of Swendsen & Wang (1987) addresses this problem by forming clusters of neighbouring pixels, then updating all of the labels within a cluster to the same value. When simulating from the prior, such as a Potts model without an external field, this algorithm is very efficient.
SW function in the PottsUtils package is implemented in a combination of
swNoData function in bayesImageS is implemented using RcppArmadillo, which gives it a speed advantage. It is worth noting that the intention of bayesImageS is not to replace PottsUtils. Rather, an efficient Swendsen-Wang algorithm is used as a building block for implementations of ABC (Grelaud et al., 2009), path sampling (Gelman & Meng, 1998), and the exchange algorithm (Murray et al., 2006). These other algorithms will be covered in future posts.
I’ve been quiet for the last few months because the work I was doing on pseudolikelihood and critical temperature developed into a paper that I had to write up and submit. Hopefully this will be the second journal article of my PhD thesis by publication, not counting the short paper I’ve submitted to the ICCR 2013 proceedings, nor the 5 manuscript pages of lit review that I’ve contributed to a middle-author paper that’s currently in prep.
A reader of this blog sent me a question about my posts on DICOM (part 1, part 2): specifically how to interact with the DICOM image as a stack of 2D slices. There are two components of this: first, I store the voxel intensity values as a numeric 3D array in R; secondly, I use the excellent misc3d package by Dai Feng and Luke Tierney to scroll interactively through the axial slices. Details (and R code) after the jump:
In my previous post I demonstrated how to read image data and meta-data from DICOM files, using the R package oro.dicom. In radiotherapy there is another kind of data that I work with, known as a DICOM-RT structure set. This is a special type of DICOM file that contains geometric contours of regions of interest, rather than image data.
The R packages oro.dicom and oro.nifti  provide the ability to import and visualise medical images that have been stored as DICOM  or ANALYZE™/NIfTI  files. This is handy because DICOM is widely used for clinical data, while NIfTI is a more compact format that stores only 2 files per 3D volume, rather than dozens.
library(oro.dicom) dcmImages <- readDICOM("/dev/DICOMRT2", verbose = TRUE, recursive = FALSE, exclude = "sql")
## 412 files to be processed by readDICOM()