gtDefFwdProjGvdm2NW.m

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 = gtAcqGetOmegaStep(parameters, det_ind);

    g = gtMathsRod2OriMat(grain.R_vector);

    if (strcmpi(fedpars.defmethod, 'rod_rightstretch'))
        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, :);
            pls_orig = gtVectorLab2Cryst(pls_orig, g);
        end
    else
        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('gtDefFwdProjGvdm2W: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 (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);

    % 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);

    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

        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);

        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; 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