forked from JeffFessler/mirt
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathfbp_fan_arc.m
232 lines (199 loc) · 6.4 KB
/
fbp_fan_arc.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
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
function img = fbp_fan_arc(sino, G, window, varargin)
%function img = fbp_fan_arc(sino, G, window, [options])
%|
%| FBP reconstruction of fan-beam tomography data collected with arc detector.
%| Assumes system 14 geometry and extracts parameters from Gtomo2_dsc object.
%| Assumes that the arc is focused on the source, e.g., 3rd generation X-ray CT.
%| See fbp_fan_arc_example.m for example.
%|
%| in:
%| sino [nb,na,nz] sinogram(s) (line integrals)
%| G [nd,np] Gtomo2_dscmex object, or empty if varargin used!
%| window [npad] 'ramp', or 'hann', or array (default: 'ramp')
%| if array, then use samples [-K/2, K/2)
%| options:
%| 'ia_skip' [int] downsample in angle to save time for quick tests
%| 'rmax' [value] maximum radius within FOV
%| (default: usual fully-sampled FOV)
%| 'mask'
%| 'nx'
%| 'ny'
%| 'orbit'
%| 'orbit_start'
%| 'pixel_size'
%| 'ds'
%| 'strip_width'
%| 'center_x'
%| 'center_y'
%| 'source_offset'
%| 'channel_offset'
%| 'flip_x'
%| 'flip_y'
%|
%| 'dsd = []; % dis_src_det
%| 'dod = []; % dis_iso_det
%| 'obj2det_y = [];
%|
%| out:
%| img [nx,ny,nz] reconstructed image(s)
%|
%| source_offset has the same units (e.g., mm) as pixel_size etc.
%| channel_offset is an integer or fraction thereof, e.g., 0.25.
%| channel offset is relative to centerline between two central channels.
%|
%| References:
%| Gullberg et al, T-MI, Mar. 1986 (gullberg:86:raf).
%| Kak & Slaney, p. 77. Fessler tomography chapter.
%|
%| Copyright December 2001. Idris Elbakri, The University of Michigan
%| Accelerated by Jeff Fessler 2002-2-19, rewritten 2004-5-21.
if nargin == 1 && streq(sino, 'test'), run_mfile_local('fbp_fan_arc_example'), return, end
if nargin < 2, ir_usage, end
if ~isvar('window') || isempty(window)
arg.window = 'ramp';
else
arg.window = window;
end
% defaults
arg.mask = [];
arg.nx = 2^9;
arg.ny = 2^9;
arg.ia_skip = 1;
arg.rmax = [];
arg.orbit = pi;
arg.orbit_start = 0;
arg.pixel_size = 1;
arg.ds = 1;
arg.strip_width = []; % not needed really
arg.center_x = 0;
arg.center_y = 0;
arg.source_offset = 0;
arg.channel_offset = 0;
arg.flip_x = 1;
arg.flip_y = 1;
arg.dsd = []; % dis_src_det
arg.dod = []; % dis_iso_det
arg.obj2det_y = [];
arg = vararg_pair(arg, varargin, 'subs', ...
{'src_det_dis', 'dsd'; 'Dsd', 'dsd';
'obj2det_x', 'dod'; 'Dod', 'dod'});
% extract arguments from system 14 Gtomo2_dsc object, if provided
if ~isempty(G)
arg = fbp_fan_arc_parse_struct(arg, G);
if size(sino,1) ~= G.arg.nb || size(sino,2) ~= G.arg.na
error 'bad sino size'
end
end
if arg.obj2det_y ~= arg.dod, error 'only circular orbit implemented', end
img = fbp_fan_arc_do(sino, arg.window, ...
arg.mask, arg.nx, arg.ny, ...
arg.ia_skip, arg.rmax, ...
arg.orbit, arg.orbit_start, ...
arg.pixel_size, arg.ds, arg.strip_width, ...
arg.center_x, arg.center_y, ...
arg.source_offset, arg.channel_offset, ...
arg.flip_x, arg.flip_y, ...
arg.dsd, arg.dod);
function img = fbp_fan_arc_do(sino, window, ...
mask, nx, ny, ...
ia_skip, rmax, ...
orbit, orbit_start, ...
pixel_size, ds, strip_width, ...
center_x, center_y, ...
source_offset, channel_offset, ...
flip_x, flip_y, ...
dsd, dod);
[nb na nz] = size(sino);
%
% fan beam FBP step 1: weight the sinogram(s)
%
nn = [-(nb-1)/2:(nb-1)/2]' - channel_offset;
dso = dsd - dod; % src to "isocenter" distance
ss = nn * ds; % sample locations
w1 = abs(dso * cos(ss/dsd) - source_offset * sin(ss/dsd)) / dsd; % 1D weighting
sino = sino .* repmat(w1, [1 na nz]);
clear nn w1 ss
%
% fan beam FBP step 2: filter the sinogram
% trick: extra zero column saves linear interpolation indexing within loop!
%
sino = fbp2_sino_filter('arc', sino, 'ds', ds, 'dsd', dsd, ...
'extra', 1, 'window', window);
%
% fan beam FBP step 3: backproject the filtered sinogram
%
% precompute as much as possible
wx = (nx+1)/2 - center_x;
wy = (ny+1)/2 - center_y;
[xc yc] = ndgrid(flip_x * ([1:nx]-wx) * pixel_size, ...
-flip_y * ([1:ny]-wy) * pixel_size);
rr = sqrt(xc.^2 + yc.^2); % [nx,ny]
smax = ((nb-1)/2-abs(channel_offset)) * ds;
if isempty(rmax)
rmax = dso * sin(smax / dsd); % alpha_max
end
mask = mask & (rr < rmax);
xc = xc(mask(:)); % [np] pixels within mask
yc = yc(mask(:));
clear wx wy rr smax
ia = [0:na-1]';
betas = orbit_start + orbit * ia / na; % [na]
wb = (nb+1)/2 + channel_offset;
img = 0;
% loop over each projection angle
for ia=1:ia_skip:na
ticker(mfilename, ia, na)
beta = betas(ia);
d_loop = dso + xc * sin(beta) - yc * cos(beta); % dso - y_beta
r_loop = xc * cos(beta) + yc * sin(beta) - source_offset; % x_beta-roff
w2 = dsd^2 ./ (d_loop.^2 + r_loop.^2); % [np] image-domain weighting
sprime_ds = (dsd/ds) * atan2(r_loop, d_loop); % s' / ds
bb = sprime_ds + wb; % [np] bin "index"
% nearest neighbor interpolation:
% ib = round(bb);
% if any(ib < 1 | ib > nb), error 'bug', end
% % trick: make out-of-sinogram indices point to those extra zeros
%% ib(ib < 1 | ib > nb) = nb+1;
% img = img + sino(ib, ia) ./ L2;
% linear interpolation:
il = floor(bb); % left bin
ir = 1+il; % right bin
% deal with truncated sinograms
ig = il >= 1 & ir <= nb;
il(~ig) = nb+1;
ir(~ig) = nb+1;
% if any(il < 1 | il >= nb), error 'bug', end
wr = bb - il; % left weight
wl = 1 - wr; % right weight
if nz > 1
wr = repmat(wr, [1 nz]);
wl = repmat(wl, [1 nz]);
img = img + (wl .* squeeze(sino(il, ia, :)) ...
+ wr .* squeeze(sino(ir, ia, :))) .* repmat(w2, [1 nz]);
else
img = img + (wl .* sino(il, ia) + wr .* sino(ir, ia)) .* w2;
end
end
img = pi / (na/ia_skip) * embed(img, mask);
% fbp_fan_arc_parse_struct()
function arg = fbp_fan_arc_parse_struct(arg, G);
arg_num = @(G, arg) str2num(arg_get(G.arg.args, arg));
arg_num_def = @(G, arg, def) str2num(arg_get(G.arg.args, arg, def));
arg.orbit = arg_num(G, 'orbit') * pi / 180; % to radians
arg.orbit_start = arg_num(G, 'orbit_start') * pi / 180; % to radians
arg.pixel_size = arg_num(G, 'pixel_size');
arg.ds = arg_num(G, 'ray_spacing'); % sample spacing [distance]
arg.strip_width = arg_num(G, 'strip_width');
arg.dsd = arg_num(G, 'src_det_dis'); % dis_src_det
arg.dod = arg_num(G, 'obj2det_x'); % dis_iso_det
arg.obj2det_y = arg_num(G, 'obj2det_y');
arg.center_x = arg_num_def(G, 'center_x', '0'); % center offset in pixel
arg.center_y = arg_num_def(G, 'center_y', '0');
arg.source_offset = arg_num(G, 'source_offset'); % r_off
arg.channel_offset = arg_num(G, 'channel_offset');
arg.flip_x = arg_num_def(G, 'flip_x', '1');
arg.flip_y = arg_num_def(G, 'flip_y', '1');
%arg.nx = G.nx;
%arg.ny = G.ny;
arg.mask = G.arg.mask;
[arg.nx arg.ny] = size(arg.mask);