ADJUST: A Dictionary-based Joint reconstruction and Unmixing method for Spectral Tomography
ADJUST is a MATLAB solver for large scale spectral tomographic inverse problems where a dictionary of spectral response of various materials is available. This algorithm jointly reconstructs and unmixes (i.e., material decomposition) spectral tomographic measurements. In this work package the following algorithms are included:
-
Unmixing-then-reconstruction (UR), Reconstruction-then-Unmixing (RU) and classic Joint reconstruction algorithm (cJoint) solves the following problem:
-
ADJUST (dictonary-based method):
The matrix W is a tomography operator of size m x n, Y are the tomographic measurements of size m x c, A contains the spatial information of materials (size: n x k), F contains the spectral information of materials (size: k x c), T is a spectral dictionary of p materials, while R is dictionary coefficient matrix. Here, m are the number of tomographic measurements, n is the size of image, c are the spectral channels, and k are the number of materials.
To run the examples on MATLAB (r2019 and above recommended), please install the following packages:
- ASTRA Toolbox: https://github.com/astra-toolbox/astra-toolbox
- SPOT operator: https://github.com/mpf/spot
- MinConf package: https://www.cs.ubc.ca/~schmidtm/Software/minConf.html
- 3D Shepp-Logan phantom package (for the 3D example only): https://www.mathworks.com/matlabcentral/fileexchange/9416-3d-shepp-logan-phantom
Eight examples are provided with the code. Before running any example, please run the startup.m script first.
- example1_SheppLogan: The classic Shepp-Logan phantom with variable size
- example2_Disk: A custom-made disk example with variable size and number of disks
- example3_Thorax: The synthetic phantom mimics internal body structure in an abstract fashion (taken from the CONRAD software framework: https://git5.cs.fau.de/PyConrad/pyCONRAD)
- example4_SheppLogan_SparseAngle: The Shepp-Logan phantom with only 10 projection angles (adjustable)
- example5_SheppLogan_LimitedView: The Shepp-Logan phantom with a missing wedge of 60 degrees (adjustable)
- example6_MixedDisk: A custom-made disk example with variable size and number of disks, where the disks consists of multiple materials
- example7_SheppLogan3D: 3D version of the Shepp-Logan phnatom with adjustable size
- example8_microCTData.m: A demonstration of the approach on real micro-CT data from a natural rock sample. The source of the data can be found further below.
On each case the four implemented algorithms (UR, RU, cJoint and ADJUST) can be tested. The perfomance of these algorithms are captured by computing and presenting the MSE, PSNR, and the SSIM between the reconstructed material spatial maps and the ground truth.
As an example, below are results for a numerical study on the Shepp-Logan phantom:
For the micro-CT dataset, four different materials (lead, tungsten, quartz and gold) can be reconstructed. An example reconstruction is shown below:
To reproduce the provided material spectra matrices (F and T), the scripts for generating these are provided in the python_spectral folder. These scripts are written in Python (version 3.7 or higher). Apart from this, physdata, csv and scipy are necessary to run these scripts.
The algorithms implemented in this MATLAB package are described in following paper. If you use (parts of) this code in a publication, we would appreciate it if you would refer to:
@article{
title={ADJUST: A Dictionary-Based Joint Reconstruction and Unmixing Method for Spectral Tomography},
author={Zeegers, Math{\'e} T and Kadu, Ajinkya and van Leeuwen, Tristan and Batenburg, Kees Joost},
journal={Inverse Problems},
year={2022},
volume={22},
number={12},
pages={125002},
publisher={IOP Publishing},
doi={10.1088/1361-6420/ac932e}
}
The preprint can be found here.
The micro-CT dataset required to run the eighth script is introduced in this paper by Sittner et al. (2020). The dataset can be found here.
Code written by:
- Ajinkya Kadu (aak [at] cwi [dot] nl)
- Mathé Zeegers (m [dot] t [dot] zeegers [at] cwi [dot] nl).