Computed Tomography using CUQIpy#
Overview#
This benchmark focuses on image reconstruction in X-ray computed tomography (CT). In CT, X-ray are passed through an object of interest and projection images are recorded at all orientations. Materials of different density absorb X-rays by different amounts quantified by the linear attenuation coefficients. A linear forward model describes how the 2D linear attenuation image of an object results in a measured set of projections, typically known as a sinogram. The benchmark uses the library CUQIpy to specify the linear forward model, based on a matrix constructed with the library AIR Tools II. It defines a posterior distribution for a 2D X-ray CT image reconstruction problem problem, with a Gaussian noise distribution ie likelihood and four different choices of 2D prior distributions with configurable parameters.
Plot of exact solution, noise-free data, noisy data and noise:
True image |
Noise-free data |
---|---|
Noisy data |
Noise |
Mean and standard deviation plots for posterior samples using different priors:
Prior |
Posterior mean |
Posterior std |
---|---|---|
Gaussian \(\delta=0.01\) |
||
GMRF \(\delta=0.01\) |
||
LMRF \(\delta=0.1\) |
||
CMRF \(\delta=0.01\) |
Run#
docker run -it -p 4243:4243 linusseelinger/benchmark-cuqi-ct
Properties#
Model |
Description |
---|---|
CT_Gaussian |
Posterior distribution with Gaussian prior |
CT_GMRF |
Posterior distribution with Gaussian Markov Random Field prior |
CT_LMRF |
Posterior distribution with Laplacian Markov Random Field prior |
CT_CMRF |
Posterior distribution with Cauchy Markov Random Field prior |
CT_ExactSolution |
Exact solution to the CT problem |
CT_Gaussian#
Mapping |
Dimensions |
Description |
---|---|---|
input |
[256**2] |
Signal \(\mathbf{x}\) |
output |
[1] |
Log PDF \(\pi(\mathbf{b}\mid\mathbf{x})\) of posterior with Gaussian prior |
Feature |
Supported |
---|---|
Evaluate |
True |
Gradient |
True |
ApplyJacobian |
False |
ApplyHessian |
False |
Config |
Type |
Default |
Description |
---|---|---|---|
delta |
double |
0.01 |
The prior parameter \(\delta\) (see below). |
CT_GMRF#
Mapping |
Dimensions |
Description |
---|---|---|
input |
[256**2] |
Signal \(\mathbf{x}\) |
output |
[1] |
Log PDF \(\pi(\mathbf{b}\mid\mathbf{x})\) of posterior with Gaussian Markov Random Field prior |
Feature |
Supported |
---|---|
Evaluate |
True |
Gradient |
True |
ApplyJacobian |
False |
ApplyHessian |
False |
Config |
Type |
Default |
Description |
---|---|---|---|
delta |
double |
0.01 |
The prior parameter \(\delta\) (see below). |
CT_LMRF#
Mapping |
Dimensions |
Description |
---|---|---|
input |
[256**2] |
Signal \(\mathbf{x}\) |
output |
[1] |
Log PDF \(\pi(\mathbf{b}\mid\mathbf{x})\) of posterior Laplacian Markov Random Field prior |
Feature |
Supported |
---|---|
Evaluate |
True |
Gradient |
False |
ApplyJacobian |
False |
ApplyHessian |
False |
Config |
Type |
Default |
Description |
---|---|---|---|
delta |
double |
0.01 |
The prior parameter \(\delta\) (see below). |
CT_CMRF#
Mapping |
Dimensions |
Description |
---|---|---|
input |
[256**2] |
Signal \(\mathbf{x}\) |
output |
[1] |
Log PDF \(\pi(\mathbf{b}\mid\mathbf{x})\) of posterior with Cauchy Markov Random Field prior |
Feature |
Supported |
---|---|
Evaluate |
True |
Gradient |
True |
ApplyJacobian |
False |
ApplyHessian |
False |
Config |
Type |
Default |
Description |
---|---|---|---|
delta |
double |
0.01 |
The prior parameter \(\delta\) (see below). |
CT_ExactSolution#
Mapping |
Dimensions |
Description |
---|---|---|
input |
[0] |
No input to be provided. |
output |
[256**2] |
Returns the exact solution \(\mathbf{x}\) for the deconvolution problem. |
Feature |
Supported |
---|---|
Evaluate |
True |
Gradient |
False |
ApplyJacobian |
False |
ApplyHessian |
False |
Config |
Type |
Default |
Description |
---|---|---|---|
None |
Mount directories#
Mount directory |
Purpose |
---|---|
None |
Source code#
The benchmark is set up in the file server.py as well as data_script.py in the same folder as this readme here. The benchmark uses CUQIpy.
Description#
The CT problem is defined by the inverse problem
where \(\mathbf{b}\) is an \(m\)-dimensional random vector representing the observed (vectorized) CT sinogram data, \(\mathbf{A}\) is an \(m\times n\) matrix representing the CT projection operator or “system matrix”, \(\mathbf{x}\) is an \(n\)-dimensional random vector representing the unknown (vectorized) square image, and \(\mathbf{e}\) is an \(m\)-dimensional random vector representing the noise.
The matrix \(\mathbf{A}\) is generated by using the MATLAB package AIR Tools II, described in the publication. Specifically, the following code was used to generate the forward model matrix:
[A, b, x] = paralleltomo(256, 0:6:174, N);
save A256_30.mat A
This corresponds to an image of \(256 \times 256\) pixels, hence \(n=256^2=65536\), scanned in a parallel-beam geometry with 30 projections taking over 180 degrees with angular steps of 6 degrees, each projection consisting of 256 X-ray lines, i.e., \(m = 30*256=7680\). See the above links for more details.
This benchmark defines a posterior distribution over \(\mathbf{x}\) given \(\mathbf{b}\) as
where \(\pi(\mathbf{b}|\mathbf{x})\) is a likelihood function and \(\pi(\mathbf{x})\) is a prior distribution.
The noise is assumed to be Gaussian with a known noise level, and so the likelihood is defined via
where \(\mathbf{I}_m\) is the \(m\times m\) identity matrix and \(\sigma=0.01\) defines the noise level.
The prior can be configured by choosing of the following assumptions about \(\mathbf{x}\), where we here consider x a square image with pixels \(x_{i,j}\) for \(i,j = 1, \dots, 256\):
Gaussian
: Gaussian (Normal) distribution: \(x_{i,j} \sim \mathcal{N}(0, \delta)\).GMRF
Gaussian Markov Random Field: \(x_{i,j}-x_{i-1,j} \sim \mathcal{N}(0, \delta)\) and \(x_{i,j}-x_{i,j-1} \sim \mathcal{N}(0, \delta)\)CMRF
Cauchy Markov Random Field: \(x_{i,j}-x_{i-1,j} \sim \mathcal{C}(0, \delta)\) and \(x_{i,j}-x_{i,j-1} \sim \mathcal{C}(0, \delta)\)LMRF
Laplace Markov Random Field: \(x_{i,j}-x_{i-1,j} \sim \mathcal{L}(0, \delta)\) and \(x_{i,j}-x_{i,j-1} \sim \mathcal{L}(0, \delta)\)
where \(\mathcal{C}\) is the Cauchy distribution and \(\mathcal{L}\) is the Laplace distribution. The parameter \(\delta\) is the prior parameter and is configurable (see above).
The choice of prior is specified by providing the name to the HTTP model. In this case CT_Gaussian
, CT_GMRF
, CT_LMRF
, and CT_CMRF
, respectively. See um-bridge Clients for more details.
In addition to the HTTP models for the posterior, there is also an HTTP model for the exact solution to the problem. This model is called CT_ExactSolution
and returns exact phantom used to generate the synthetic data when called.
Using CUQIpy the CT benchmark is defined using the file data_script.py and server.py provided here.