Newer
Older
Nicola Vigano
committed
function shape_funcs = gtDefShapeFunctionsFwdProjUVW(sampler, factor, dtype, recenter_sf)
Nicola Vigano
committed
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
if (~exist('dtype', 'var') || isempty(dtype))
dtype = 'single';
end
if (~exist('recenter_sf', 'var') || isempty(recenter_sf))
recenter_sf = true;
end
interpolated_voxel = false;
det_ind = sampler.detector_index;
labgeo = sampler.parameters.labgeo;
% samgeo = sampler.parameters.samgeo;
% recgeo = sampler.parameters.recgeo(det_ind);
% detgeo = sampler.parameters.detgeo(det_ind);
labgeo.rotcomp = gtMathsRotationMatrixComp(labgeo.rotdir', 'col');
sampler.parameters.labgeo = labgeo;
fed_pars_detector = struct(...
'blobsizeadd', [0 0 0], 'psf', [], 'apply_uv_shift', recenter_sf);
fedpars = struct( ...
'dcomps', [1 1 1 0 0 0 0 0 0] == 1, ...
'defmethod', 'rod_rightstretch', ...
'detector', repmat( fed_pars_detector, [numel(sampler.parameters.detgeo), 1]));
voxel_factor = factor([1 1 1]);
num_sub_voxels = prod(voxel_factor);
ones_nsv = ones(num_sub_voxels, 1);
gv = struct('cs', [], 'ind', [], 'dm', [], 'pow', []);
[gv.ind(:, 1), gv.ind(:, 2), gv.ind(:, 3)] = ind2sub(voxel_factor, 1:num_sub_voxels);
gv.ind = gv.ind';
ref_gr = sampler.get_reference_grain();
gv.cs = ref_gr.center(ones_nsv, :)';
% computing the subsampling, used to deterine the shape functions
voxel_size = sampler.stats.sampling.gaps;
if (interpolated_voxel)
half_voxel_size = voxel_size; % Enlarging to neighbours
steps = voxel_size ./ voxel_factor * 2;
else
half_voxel_size = voxel_size ./ 2;
steps = voxel_size ./ voxel_factor;
end
half_steps = steps ./ 2;
beg_pos = - half_voxel_size + half_steps;
end_pos = + half_voxel_size - half_steps;
x_steps = beg_pos(1):steps(1):end_pos(1);
y_steps = beg_pos(2):steps(2):end_pos(2);
z_steps = beg_pos(3):steps(3):end_pos(3);
num_steps = [numel(x_steps), numel(y_steps), numel(z_steps)];
ones_x_steps = ones(1, num_steps(1));
gv.d = zeros(3, prod(num_steps));
pos = 0;
for ii_z = z_steps
for ii_y = y_steps
gv.d(:, (pos+1):(pos+num_steps(1))) ...
= [ x_steps; ones_x_steps * ii_y; ones_x_steps * ii_z];
pos = pos + num_steps(1);
end
end
if (interpolated_voxel)
gv.pow = prod(1 - abs(gv.d ./ half_voxel_size(ones_nsv, :)'), 1);
else
gv.pow = ones(1, num_sub_voxels, dtype);
end
gv.pow = gv.pow / sum(gv.pow);
gv.used_ind = find(gv.pow > 0);
nuv = numel(gv.used_ind);
ones_nuv = ones(1, nuv);
sel_bls = find(sampler.selected);
inc_bls = sampler.ondet(sampler.included(sampler.selected));
% computing the size of the blobs
shape_funcs_size(:, 3) = ceil(2 * atand(sampler.stats.sampling.gaps(3)) ...
/ labgeo.omstep .* ref_gr.allblobs.lorentzfac(inc_bls)) + 1;
shape_funcs_size(:, 1) = sampler.stats.proj.max_pixel_dist_u;
shape_funcs_size(:, 2) = sampler.stats.proj.max_pixel_dist_v;
if (interpolated_voxel)
shape_funcs_size = shape_funcs_size * 2;
end
shape_funcs_size = 2 * ceil(shape_funcs_size) + 1; % Imposing odd shape! -> for centering
half_sf_size = floor(shape_funcs_size / 2);
% Only one lattice allowed for the moment!
num_ors = numel(sampler.lattice(1).gr);
shape_funcs = cell(num_ors, 1);
nbl = numel(sel_bls);
fprintf('Computing UVW shape functions: ')
c = tic();
for ii_g = 1:num_ors
num_chars = fprintf('%03d/%03d', ii_g, num_ors);
or = sampler.lattice(1).gr{ii_g};
bl = repmat(gtFedBlobDefinition(), nbl, 1);
or_uvw = or.allblobs.detector(det_ind).uvw(sel_bls, :);
for ii_b = 1:nbl
% Load diffraction paramaters of blobs
% Initial strained non-signed plane normal in SAMPLE reference
bl(ii_b).plorig = or.allblobs.plorig(sel_bls(ii_b), :);
bl(ii_b).sinth0 = or.allblobs.sintheta(sel_bls(ii_b)); % initial sin(theta)
bl(ii_b).omega0 = or.allblobs.omega(sel_bls(ii_b)); % initial omega
bl(ii_b).omind = or.allblobs.omind(sel_bls(ii_b)); % initial omega index (1..4)
% Initial centroid in image
bl(ii_b).u0im = or_uvw(ii_b, :);
% Projection information!!
bl(ii_b).bbuim = round(or_uvw(ii_b, 1) + [-half_sf_size(ii_b, 1), half_sf_size(ii_b, 1)]);
bl(ii_b).bbvim = round(or_uvw(ii_b, 2) + [-half_sf_size(ii_b, 2), half_sf_size(ii_b, 2)]);
bl(ii_b).bbwim = round(or_uvw(ii_b, 3) + [-half_sf_size(ii_b, 3), half_sf_size(ii_b, 3)]);
bl(ii_b).bbsize = [ ...
(bl(ii_b).bbuim(2) - bl(ii_b).bbuim(1) + 1), ...
(bl(ii_b).bbvim(2) - bl(ii_b).bbvim(1) + 1), ...
(bl(ii_b).bbwim(2) - bl(ii_b).bbwim(1) + 1) ];
% Blob origin and center [u, v, w]:
bl(ii_b).bbor = [bl(ii_b).bbuim(1)-1, bl(ii_b).bbvim(1)-1, bl(ii_b).bbwim(1)-1];
bl(ii_b).bbcent = [ ...
round( (bl(ii_b).bbuim(1) + bl(ii_b).bbuim(2) ) / 2 ), ...
round( (bl(ii_b).bbvim(1) + bl(ii_b).bbvim(2) ) / 2 ), ...
round( (bl(ii_b).bbwim(1) + bl(ii_b).bbwim(2) ) / 2 ) ];
% u in blob = u in image minus u at blob origin
bl(ii_b).u0bl = bl(ii_b).u0im - bl(ii_b).bbor;
bl(ii_b).uv_shift = bl(ii_b).u0im(1:2) - bl(ii_b).bbcent(1:2);
bl(ii_b).int = zeros(bl(ii_b).bbsize, dtype);
bl(ii_b).intm = [];
end
for ii_b = nbl:-1:1
gvb(ii_b).uc = zeros(3, nuv, dtype);
gvb(ii_b).ucd = zeros(3, nuv, dtype);
gvb(ii_b).ucblmod = zeros(3, nuv, dtype);
gvb(ii_b).ucbl8 = zeros(1, nuv*8, dtype);
gvb(ii_b).ucbl8ini= zeros(1, nuv*8, dtype);
gvb(ii_b).ucbl8ok = 0;
% Initial u of gv centre relative to blob origin
gvb(ii_b).uc0bl = cast(bl(ii_b).u0im - bl(ii_b).bbor, dtype)';
gvb(ii_b).uc0bl = gvb(ii_b).uc0bl(:, ones_nuv);
end
[gvb, bl] = gtFedFwdProjExact(gv, gvb, bl, fedpars, sampler.parameters, det_ind, false);
% Cutting borders
for ii_b = 1:nbl
diff_inds_w = [...
find(sum(sum(bl(ii_b).int, 1), 2), 1, 'first') - 1, ...
find(sum(sum(bl(ii_b).int, 1), 2), 1, 'last') - bl(ii_b).bbsize(3); ...
];
bl(ii_b).bbwim = bl(ii_b).bbwim + diff_inds_w;
bl(ii_b).mbbu = bl(ii_b).bbuim;
bl(ii_b).mbbv = bl(ii_b).bbvim;
bl(ii_b).mbbw = bl(ii_b).bbwim;
inds_w = [1, bl(ii_b).bbsize(3)] + diff_inds_w;
bl(ii_b).intm = bl(ii_b).int(:, :, inds_w(1):inds_w(2));
bl(ii_b).int = [];
bl(ii_b).bbsize(3) = size(bl(ii_b).intm, 3);
end
shape_funcs{ii_g} = bl;
fprintf(repmat('\b', [1, num_chars]));
end
fprintf('Done in %g seconds.\n', toc(c))
end