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'))
pls_orig = grain.allblobs(det_ind).plcry(ref_sel, :);
else
% We bring the plane normals back to the status of pl_cry
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
pls_orig = grain.allblobs(det_ind).plorig(ref_sel, :);
end
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
sinths = grain.allblobs(det_ind).sintheta(ref_sel);
ominds = grain.allblobs(det_ind).omind(ref_sel);
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
if (isfield(fedpars, 'ref_omega') && ~isempty(fedpars.ref_omega))
ref_oms = fedpars.ref_omega(ref_sel);
else
ref_oms = grain.allblobs(det_ind).omega(ref_sel);
end
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
ref_ns = grain.allblobs(det_ind).eta(ref_sel);
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
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];
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