-
Notifications
You must be signed in to change notification settings - Fork 7
/
Copy pathdemo_polyquant_fanbeam.m
74 lines (71 loc) · 3.59 KB
/
demo_polyquant_fanbeam.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Fanbeam Polyquant demo
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Description
% ----------
% This script demonstrates direct quantitative reconstruction into relative
% electron density from polyenergetic X-ray CT measurements. There are 6
% samples, including various anatomical sites from the ICRP 110 female
% computational phantom, and one ('implant') with added double titanium
% hips.
% The measurements were generated using Gate software and consist of:
% proj -- the total measured photons after detector response function.
% scat -- the photons from scatter (very small and ignored here).
% i0 -- the incident source flux.
% The structure in fandata.specData contains the spectral information of
% the source (only a subsampling of source used in simulation):
% specData.energy -- the energies (MeV) in the subsampled spectrum.
% specData.spectrum -- the subsampled source spectrum.
% specData.response -- the detector response function.
% specData.hinge -- the location of the piecewise linear fit
% transitions, for 3 linear sections.
% specData.knee -- contains the equations for the piecewise linear fits
% between relative electron density and each energy in
% specData.energy. This was fitted against the
% biological materials in the ICRP 89 and for titanium
% (density = 4.506 g/cm3).
% Things to try
% ------------
% o Test the various samples by changing 'sample' string.
% o Test the influence on regularisation constant 'lambda'.
% o Try removing the 'mode.proxFun' line, to give ML reconstruction.
% o Adjust 'mode.maxIter', 'mode.tau' and 'mode.nSplit' to change
% convergence properties.
% o Set 'mode.numLinFit = 1' to see beam-hardening artefacts, which is
% equivalent to a linearised model (monoenergetic).
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Created: 26/04/2019
% Last edit: 31/05/2019
% Jonathan Hugh Mason
%
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% References: (please cite if making use of this code or its methods)
% Jonathan H Mason et al 2017 Phys. Med. Biol. 62 8739
% Jonathan H Mason et al 2018 Phys. Med. Biol. 63 225001
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%% Load data
load data/attenuationDb
load data/fandata
addpath(genpath('.'));
%% Setup the geometry
ig = image_geom('nx',137,'ny',299,'dx',0.1775,'dy',0.1775);
sg = sino_geom('fan','dfs',inf,'dsd',120.3,'dso',64.5,'units','cm',...
'ns',512,'ds',85/512,'strip_width','ds','na',360,'down',1);
A = Gtomo2_dscmex(sg, ig);
%% Polyquant setup
sample = 'chest'; % 'brain', 'head', 'abdomen', 'pelvis', 'implant'
mode = [];
mode.useConst = true; % just to offset objective function to better range
mode.contrast = [0.5,1.4]; % display contrast
mode.verbose = 2; % 0 = no output; 1 = text output; 2 = text+image output
mode.tau = 5; % conservative choice
mode.nSplit = 24;
mode.maxIter = 500; % more iterations recommended for implant
lambda = 0.5; % can be optimised for better results
mode.numLinFit = 2; % can be set to 2 for tissue, but use 3 for implant
mode.proxFun = @(z,t) prox_tv_nn(z,t*lambda);
mode.regFun = @(z) norm_tv(z);
xTrue = mat_to_den(attenuationDb,single(fandata.(sample).mat));
%% Perform the reconstruction
out = polyquant(mode,fandata.specData,fandata.(sample).proj,fandata.i0,A,xTrue);
fprintf('Reconstructed with PSNR = %.2f dB\n',20*log10(max(xTrue(:))./out.rmse(end)));