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