/
FDKReconstruction.m
150 lines (121 loc) · 5.26 KB
/
FDKReconstruction.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
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This matlab script computes an FDK reconstruction for one of the data sets
% described in
% "A Cone-Beam X-Ray CT Data Collection Designed for Machine Learning" by
% Henri Der Sarkissian, Felix Lucka, Maureen van Eijnatten,
% Giulia Colacicco, Sophia Bethany Coban, Kees Joost Batenburg
% Please note that the reconstruction shown in the paper were performed
% using the corresponding Python code, not this matlab code!
%
% author: Felix Lucka
% date: 18.03.2019
% last update: 06.11.2020
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
clc
close all
% add ASTRA toolbox to the matlab path
% Note that for obtaining a comparable scaling of the image intensities
% between FDK and iterative reconstructions, it is required to use a
% development version of the ASTRA toolbox more recent than 1.9.0dev!
addpath(genpath('~/astra-toolbox/matlab'))
% the script assumes that all zip files were extracted into one folder,
% please enter the path on your system here
data_dir_root = '~/Walnuts/';
%% set reconstruction parameters
% which walnut data should be reconstructed?
walnut_id = 1;
% define which source orbit data should be used (2 = middle orbit)
orbit_id = 2;
% define a sub-sampling factor in angular direction
angluar_sub_sampling = 1;
% number of voxels per mm in one direction (higher = larger res)
voxel_per_mm = 10;
%% set up scanning and volume geometry
% construct path to data
data_dir = [data_dir_root 'Walnut' int2str(walnut_id) ...
'/Projections/tubeV' int2str(orbit_id) '/'];
% generate reconstruction geometry
vol_sz = (50 * voxel_per_mm + 1) * ones(1, 3);
vol_geom = astra_create_vol_geom(vol_sz);
% By default, ASTRA assumes a voxel size of 1, we need to scale the
% reconstruction space here by the actual voxel size
vol_geom.option.WindowMaxX = vol_geom.option.WindowMaxX / voxel_per_mm;
vol_geom.option.WindowMinX = vol_geom.option.WindowMinX / voxel_per_mm;
vol_geom.option.WindowMaxY = vol_geom.option.WindowMaxY / voxel_per_mm;
vol_geom.option.WindowMinY = vol_geom.option.WindowMinY / voxel_per_mm;
vol_geom.option.WindowMaxZ = vol_geom.option.WindowMaxZ / voxel_per_mm;
vol_geom.option.WindowMinZ = vol_geom.option.WindowMinZ / voxel_per_mm;
% set up projection geometry
detector_sz = [972, 768]; % detector size
proj_geom = [];
proj_geom.type = 'cone_vec';
proj_geom.DetectorRowCount = detector_sz(1);
proj_geom.DetectorColCount = detector_sz(2);
proj_geom.Vectors = importdata([data_dir 'scan_geom_corrected.geom']);
% sub-sample in angle, note that the total number of projection is in fact 1201, but the
% first and last projection come from the same angle and are omitted here
proj_geom.Vectors = proj_geom.Vectors(1:angluar_sub_sampling:1200, :);
n_pro = size(proj_geom.Vectors, 1);
%% read in and normalize all data
% get all projections
pro_files = dir([data_dir, 'scan_*.tif']);
% we need to read in the projection in reverse order due to the portrait
% mode acquision
pro_files = pro_files(1201:-angluar_sub_sampling:2);
% transformation to apply to each image, we need to get the image from
% the way the scanner reads it out into to way described in the projection
% geometry
trafo = @(x) flipud(x)';
% read in dark and flat field
dark_field = trafo(double(imread([data_dir 'di000000.tif'])));
flat_field1 = trafo(double(imread([data_dir 'io000000.tif'])));
flat_field2 = trafo(double(imread([data_dir 'io000001.tif'])));
% we simply average the flat fields here
flat_field = (flat_field1 + flat_field2) / 2;
data = zeros([n_pro, size(dark_field)]);
% loop over projection data
for i_pro = 1:n_pro
pro = trafo(double(imread([data_dir, pro_files(i_pro).name])));
% flat and dark field correction
data(i_pro,:, :) = (pro - dark_field)./ (flat_field - dark_field);
end
% reset values smaller or equal to 0
data(data <= 0) = min(data(data > 0));
% values larger than 1 are clipped to 1
data = min(data, 1);
% log data
data = - log(data);
% permute data to ASTRA convention
data = permute(data, [3,1,2]);
%% compute the FDK reconstruction
% Create astra objects for the reconstruction
rec_id = astra_mex_data3d('create', '-vol', vol_geom, 0);
proj_id = astra_mex_data3d('create', '-proj3d', proj_geom, data);
clear data
% Set up the parameters for a reconstruction algorithm using the GPU
cfg = astra_struct('FDK_CUDA');
cfg.ReconstructionDataId = rec_id;
cfg.ProjectionDataId = proj_id;
% Create the algorithm object from the configuration structure
alg_id = astra_mex_algorithm('create', cfg);
% Run FDK
clock_cmp = tic;
astra_mex_algorithm('run', alg_id);
toc(clock_cmp)
% Get the result
rec = astra_mex_data3d('get', rec_id);
% Clean up. Note that GPU memory is tied up in the algorithm object,
% and main RAM in the data objects.
astra_mex_algorithm('delete', alg_id);
astra_mex_data3d('delete', rec_id, proj_id);
%% simple slice visulization
figure();
sliceX = squeeze(rec(ceil(vol_sz(1)/2), :, :));
sliceY = squeeze(rec(:, ceil(vol_sz(1)/2), :));
sliceZ = squeeze(rec(:, :, ceil(vol_sz(1)/2)));
clim = [0, max([sliceX(:);sliceY(:);sliceZ(:)])];
subplot(1, 3, 1); imagesc(sliceX, clim);
subplot(1, 3, 2); imagesc(sliceY, clim);
subplot(1, 3, 3); imagesc(sliceZ, clim);