Newer
Older
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)
om_step = gtAcqGetOmegaStep(parameters, det_ind);
nbl = numel(ref_sel);
uinds = gv.used_ind;
nv = numel(uinds);
use_diffractometer = (isfield(parameters, 'diffractometer') ...
&& numel(parameters.diffractometer) >= det_ind ...
&& det_ind > 1);
labgeo = parameters.labgeo;
samgeo = parameters.samgeo;
if (use_diffractometer)
diff_acq = parameters.diffractometer(det_ind);
diff_ref = parameters.diffractometer(1);
rotcomp = gtMathsRotationMatrixComp(diff_acq.axes_rotation', 'col');
rotdir = diff_acq.axes_rotation';
instr_t = gtGeoDiffractometerTensor(diff_acq, 'sam2lab', ...
'reference_diffractometer', diff_ref);
else
rotcomp = gtMathsRotationMatrixComp(labgeo.rotdir', 'col');
rotdir = labgeo.rotdir';
end
if (strcmpi(fedpars.defmethod, 'rod_rightstretch'))
% In this case, the deformation is relative to [0, 0, 0], so we
% need crystal plane normals
if (isfield(grain.allblobs, 'plcry'))
Nicola Vigano
committed
pls_orig = grain.allblobs(det_ind).plcry(ref_sel, :);
else
% We bring the plane normals back to the status of pl_cry
Nicola Vigano
committed
pls_orig = grain.allblobs(det_ind).plorig(ref_sel, :);
g = gtMathsRod2OriMat(grain.R_vector);
pls_orig = gtVectorLab2Cryst(pls_orig, g);
end
else
% Here the deformation is relative to the average oriantion, so we
% need the undeformed sample plane normals
Nicola Vigano
committed
pls_orig = grain.allblobs(det_ind).plorig(ref_sel, :);
if (isfield(fedpars, 'detector') ...
&& isfield(fedpars.detector(det_ind), 'blobs_w_interp') ...
&& ~isempty(fedpars.detector(det_ind).blobs_w_interp))
blobs_w_interp = fedpars.detector(det_ind).blobs_w_interp;
if (numel(blobs_w_interp) ~= nbl)
error([mfilename ':wrong_argument'], ...
'Number of "blobs_w_interp" (%d) doesn''t match with the number of blobs (%d)', ...
numel(blobs_w_interp), nbl)
end
else
blobs_w_interp = ones(nbl, 1);
end
Nicola Vigano
committed
sinths = grain.allblobs(det_ind).sintheta(ref_sel);
ominds = grain.allblobs(det_ind).omind(ref_sel);
Nicola Vigano
committed
if (~use_diffractometer)
% 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);
end
Nicola Vigano
committed
if (isfield(fedpars, 'ref_omega') && ~isempty(fedpars.ref_omega))
ref_oms = fedpars.ref_omega(ref_sel);
else
Nicola Vigano
committed
ref_oms = grain.allblobs(det_ind).omega(ref_sel);
Nicola Vigano
committed
end
Nicola Vigano
committed
w_shifts = round(ref_oms ./ om_step ./ blobs_w_interp);
if (isfield(fedpars, 'ref_eta') && ~isempty(fedpars.ref_eta))
ref_ns = fedpars.ref_eta(ref_sel);
else
Nicola Vigano
committed
ref_ns = grain.allblobs(det_ind).eta(ref_sel);
Nicola Vigano
committed
end
n_shifts = round(ref_ns ./ eta_step);
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)
% !!! use plcry and not plsam to keep omega order below !!!
pl_orig = pls_orig(ii_b, :)';
[pl_samd, drel] = gtStrainPlaneNormals(pl_orig, defT); % unit column vectors
if (use_diffractometer)
% The plane norma ls 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)
pl_samd = gtGeoSam2Lab(pl_samd', instr_t, labgeo, samgeo, true)';
end
% 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', rotdir, 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
Nicola Vigano
committed
om = gtGrainAnglesTabularFix360deg(om, ref_oms(ii_b));
w_bl = om ./ om_step;
w_bl = w_bl / blobs_w_interp(ii_b);
w_bl = w_bl - w_shifts(ii_b);
eta = gtGeoEtaFromDiffVec(pl_lab', parameters.labgeo)';
eta = gtGrainAnglesTabularFix360deg(eta, ref_ns(ii_b));
n_bl = eta ./ eta_step;
n_bl = n_bl - n_shifts(ii_b);
nw{ii_b} = [n_bl; w_bl];
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
233
234
235
236
237
238
239
240
241
242
243
244
245
246
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] * blobs_w_interp(ii_b);
bl(ii_b).interpw = blobs_w_interp(ii_b);
o.fprintf(repmat('\b', [1, num_chars]));
end
o.fprintf('Done in %g s\n', toc(t));
end