Newer
Older
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
function bl = gtDefFwdProjGvdm2NW(grain, ref_sel, gv, fedpars, parameters, eta_step, det_ind, verbose)
if (~exist('det_ind', 'var'))
det_ind = 1;
end
if (~exist('verbose', 'var'))
verbose = true;
end
o = GtConditionalOutput(verbose);
o.fprintf('Forward projection (%s):\n', mfilename)
nbl = numel(ref_sel);
uinds = gv.used_ind;
nv = numel(uinds);
labgeo = parameters.labgeo;
samgeo = parameters.samgeo;
om_step = gtGetOmegaStepDeg(parameters, det_ind);
g = gtMathsRod2OriMat(grain.R_vector);
if (strcmpi(fedpars.defmethod, 'rod_rightstretch'))
if (isfield(grain.allblobs, 'plcry'))
pls_orig = grain.allblobs.plcry(ref_sel, :);
else
% We bring the plane normals back to the status of pl_cry
pls_orig = grain.allblobs.plorig(ref_sel, :);
pls_orig = gtVectorLab2Cryst(pls_orig, g);
end
else
pls_orig = grain.allblobs.plorig(ref_sel, :);
end
sinths = grain.allblobs.sintheta(ref_sel);
ominds = grain.allblobs.omind(ref_sel);
ref_ws = grain.allblobs.detector(det_ind).uvw(ref_sel, 3);
w_shifts = round(ref_ws);
ref_ns = grain.allblobs.eta(ref_sel);
n_shifts = round(ref_ns ./ eta_step);
% The plane normals need to be brought in the Lab reference where the
% beam direction and rotation axis are defined.
% Use the Sample -> Lab orientation transformation assuming omega=0;
% (vector length preserved for free vectors)
pls_orig = gtGeoSam2Lab(pls_orig, eye(3), labgeo, samgeo, true, false);
linear_interp = ~(isfield(fedpars, 'projector') ...
&& strcmpi(fedpars.projector, 'nearest'));
gvpow = gv.pow(1, uinds);
gvd = gv.d(:, uinds);
% Deformation tensor (relative to reference state) from its components
defT = gtFedDefTensorFromComps(gvd, fedpars.dcomps, fedpars.defmethod, 0);
%%% Computation of indices
o.fprintf(' * Computing indices and bbsizes: ')
t = tic();
nw = cell(nbl, 1);
nw_min = zeros(nbl, 2);
nw_max = zeros(nbl, 2);
valid_voxels = false(nv, nbl);
for ii_b = 1:nbl
num_chars = o.fprintf('%03d/%03d', ii_b, nbl);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate new detector coordinates
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% New deformed plane normals and relative elongations (relative to
% reference state)
pl_orig = pls_orig(ii_b, :)'; % ! use plcry and not plsam to keep omega order below!
[pl_samd, drel] = gtStrainPlaneNormals(pl_orig, defT); % unit column vectors
% New sin(theta)
sinth = sinths(ii_b) ./ drel;
% Predict omega angles: 4 for each plane normal
[om, pl_lab] = gtFedPredictOmegaMultiple(pl_samd, sinth, labgeo.beamdir', ...
labgeo.rotdir', labgeo.rotcomp, ominds(ii_b));
valid_voxels(:, ii_b) = ~isnan(om);
% Delete those where no reflection occurs
if (any(isnan(om)))
inds_bad = find(isnan(om));
gvd(:, inds_bad(1:min(10, numel(inds_bad))))
warning('gtFedFwdProjExact:bad_R_vectors', ...
'No diffraction from some elements after deformation (%d over %d) for blob %d.', ...
numel(inds_bad), numel(om), ii_b)
end
nw_bl = om ./ om_step;
nw_bl = gtGrainAnglesTabularFix360deg(nw_bl, ref_ws(ii_b), parameters);
nw_bl = nw_bl - w_shifts(ii_b);
n_bl = gtGeoEtaFromDiffVec(pl_lab', parameters.labgeo)';
n_bl = gtGrainAnglesTabularFix360deg(n_bl, ref_ns(ii_b));
n_bl = n_bl ./ eta_step - n_shifts(ii_b);
nw{ii_b} = [n_bl; nw_bl];
nw_min(ii_b, :) = min(nw{ii_b}, [], 2);
nw_max(ii_b, :) = max(nw{ii_b}, [], 2);
o.fprintf(repmat('\b', [1, num_chars]));
end
o.fprintf('Done in %g s\n', toc(t));
o.fprintf(' * Computing max BBox size and feasibilities:\n')
if (linear_interp)
blob_nw_sizes = ceil(nw_max) - floor(nw_min) + 1;
blob_orig_nw_shifts = 1 - floor(nw_min);
else
blob_nw_sizes = round(nw_max) - round(nw_min) + 1;
blob_orig_nw_shifts = 1 - round(nw_min);
end
o.fprintf(' Computed Blob N sizes: [%s]\n', sprintf(' %d', blob_nw_sizes(:, 1)));
o.fprintf(' Computed Blob W sizes: [%s]\n', sprintf(' %d', blob_nw_sizes(:, 2)));
bl = gtFwdSimBlobDefinition('sf_nw', nbl);
%%% Blob projection
o.fprintf(' * Projecting volumes: ')
t = tic();
for ii_b = 1:nbl
num_chars = o.fprintf('%03d/%03d', ii_b, nbl);
% Detector coordinates U,V in blob
nw_bl = nw{ii_b} + blob_orig_nw_shifts(ii_b * ones(1, nv), :)';
% Let's now filter valid voxels
nw_bl = nw_bl(:, valid_voxels(:, ii_b))';
gvpow_v = reshape(gvpow(valid_voxels(:, ii_b)), [], 1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Sum intensities
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Distribute value over 8 pixels
if (linear_interp)
[inds_nw, ints_nw] = gtMathsGetInterpolationIndices(nw_bl, ...
gvpow_v, fedpars.bltype);
else
inds_nw = round(nw_bl);
ints_nw = gvpow_v;
end
% Accumulate all intensities in the blob voxel-wise
% 'uvw' needs to be nx1; 'ints' is now 1xnx8
try
bl(ii_b).intm = accumarray(inds_nw, ints_nw, blob_nw_sizes(ii_b, :));
catch mexc
fprintf(['\n\nERROR\nThe error is probably caused by the' ...
' projected intensities falling outside the blob volume.' ...
'\nTry increasing the blob volume padding:' ...
' fedpars.detector(det_ind).blobsizeadd, or FIXING' ...
' YOUR CODE!!\n\n'])
disp('Blob ID:')
disp(ii_b)
disp('Blob size:')
disp(blob_nw_sizes(ii_b))
disp('Min projected NW coordinate:')
fprintf('\t%g\t%g \t(rounded: %d %d )\n', ...
min(nw_bl, [], 1), min(inds_nw, [], 1))
disp('Max projected W coordinate:')
fprintf('\t%g \t(rounded: %d )\n', ...
max(nw_bl, [], 1), max(inds_nw, [], 1))
new_exc = GtFedExceptionFwdProj(mexc, det_ind, [0 0 0]);
throw(new_exc)
end
bl(ii_b).bbsize = blob_nw_sizes(ii_b, :);
bl(ii_b).bbnim = ([nw_min(ii_b, 1), nw_max(ii_b, 1)] + n_shifts(ii_b)) .* eta_step;
im_w_low_lim = w_shifts(ii_b) - blob_orig_nw_shifts(ii_b, 2) + 1;
bl(ii_b).bbwim = [im_w_low_lim, im_w_low_lim + bl(ii_b).bbsize(2) - 1];
o.fprintf(repmat('\b', [1, num_chars]));
end
o.fprintf('Done in %g s\n', toc(t));
end