gtDefFwdProjGvdm2P.m

function bl = gtDefFwdProjGvdm2P(grain, selectedph, gv, fedpars, parameters, 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)

    ph_step = gtAcqGetPhiStep(parameters, det_ind);
    nbl = gtAcqTotNumberOfImages(parameters, det_ind);

    uinds = gv.used_ind;
    nv = numel(uinds);

    labgeo = parameters.labgeo;
    samgeo = parameters.samgeo;

    % we pick the first because they are all the same!
    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'))
            pl_orig = grain.allblobs(det_ind).plcry(1, :)';
        else
            % We bring the plane normals back to the status of pl_cry
            pl_orig = grain.allblobs(det_ind).plorig(1, :);

            g = gtMathsRod2OriMat(grain.R_vector);
            pl_orig = gtVectorLab2Cryst(pl_orig, g)';
        end
    else
        % Here the deformation is relative to the average oriantion, so we
        % need the undeformed sample plane normals
        pl_orig = grain.allblobs(det_ind).plorig(1, :)';
    end

    if (~exist('selectedph', 'var') || isempty(selectedph))
        lims_basetilt = parameters.acq(det_ind).range_basetilt;
        within_lims = grain.allblobs(det_ind).phi > lims_basetilt(1) & grain.allblobs(det_ind).phi < lims_basetilt(2);
        phinds = grain.allblobs(det_ind).phind(within_lims);
        phinds_counts = histcounts(phinds, (1:5)-0.5);
        [~, phind] = max(phinds_counts);

        o.fprintf(' - selected phind: %d (counts: [%s])\n', phind, sprintf(' %d', phinds_counts))

        selectedph = grain.allblobs(det_ind).phind == phind;
    end

    sinths = grain.allblobs(det_ind).sintheta(selectedph);
    or_pw = grain.allblobs(det_ind).detector.uvpw(selectedph, 3);
    omega = grain.allblobs(det_ind).omega(selectedph);
    phinds = grain.allblobs(det_ind).phind(selectedph);

    diff_acq = parameters.diffractometer(det_ind);
    diff_ref = parameters.diffractometer(1);

    rotcomp = gtMathsRotationMatrixComp(diff_acq.axes_basetilt', 'col');
    rotdir = diff_acq.axes_basetilt';

    instr_t = gtGeoDiffractometerTensor(diff_acq, 'sam2lab', ...
        'reference_diffractometer', diff_ref, 'angles_rotation', omega);

    pw_shifts = round(or_pw)';

    linear_interp = ~(isfield(fedpars, 'projector') ...
        && isstruct(fedpars.projector) && isfield(fedpars.projector, 'interp') ...
        && strcmpi(fedpars.projector.interp, '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();

    p = cell(nbl, 1);
    p_min = zeros(nbl, 1);
    p_max = zeros(nbl, 1);

    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_samd, drel] = gtStrainPlaneNormals(pl_orig, defT); % unit column vectors

        % 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)
        pl_samd = gtGeoSam2Lab(pl_samd', instr_t(:, :, ii_b), labgeo, samgeo, true)';

        % New sin(theta)
        sinth = sinths(ii_b) ./ drel;

        % Predict omega angles: 4 for each plane normal
        ph = gtFedPredictOmegaMultiple(pl_samd, sinth, ...
            labgeo.beamdir', rotdir, rotcomp, phinds(ii_b));
        ph = mod(ph + 180, 360) - 180;

        valid_voxels(:, ii_b) = ~isnan(ph);

        % Delete those where no reflection occurs
        if (any(isnan(ph)))
            inds_bad = find(isnan(ph));
            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(ph), ii_b)
        end

        p_bl = ph ./ ph_step;
        % Transforming into offsets from average orientation
        p_bl = p_bl - pw_shifts(ii_b);

        p{ii_b} = p_bl;

        p_min(ii_b) = min(p{ii_b});
        p_max(ii_b) = max(p{ii_b});

        if (any(p_min(ii_b) > 0) || any(p_max(ii_b) < 0))
            warning([mfilename ':wrong_result'], ...
                '\nThe average orientation seems to project outside the blob!\n')
        end

        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_p_sizes = ceil(p_max) - floor(p_min) + 1;
        blob_orig_p_shifts = 1 - floor(p_min);
    else
        blob_p_sizes = round(p_max) - round(p_min) + 1;
        blob_orig_p_shifts = 1 - round(p_min);
    end

    o.fprintf('     Computed Blob Phi sizes: [%s]\n', sprintf(' %d', blob_p_sizes));

    bl = gtFwdSimBlobDefinition('sf_p', 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
        p_bl = p{ii_b}' + blob_orig_p_shifts(ii_b);

        % Let's now filter valid voxels
        p_bl = p_bl(valid_voxels(:, ii_b));
        gvpow_v = reshape(gvpow(valid_voxels(:, ii_b)), [], 1);

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Sum intensities
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Distribute value over 8 pixels

        if (linear_interp)
            [ucbl8, ints8] = gtMathsGetInterpolationIndices(p_bl, gvpow_v, fedpars.bltype);
        else
            ucbl8 = round(p_bl);
            ints8 = gvpow_v;
        end

        % Accumulate all intensities in the blob voxel-wise
        %   'uvw' needs to be nx3; 'ints' is now 1xnx8
        try
            bl(ii_b).intm = accumarray(ucbl8, ints8, [blob_p_sizes(ii_b), 1]);
        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_p_sizes(ii_b))
            disp('Min projected P coordinates:')
            fprintf('\t%g \t(rounded: %d )\n', ...
                min(p_bl, [], 1), min(ucbl8, [], 1))
            disp('Max projected P coordinates:')
            fprintf('\t%g \t(rounded: %d )\n', ...
                max(p_bl, [], 1), max(ucbl8, [], 1))
            new_exc = GtFedExceptionFwdProj(mexc, det_ind, [0 0 0]);
            throw(new_exc)
        end

        bl(ii_b).bbsize = blob_p_sizes(ii_b);

        im_low_lims = pw_shifts(ii_b) - blob_orig_p_shifts(ii_b) + 1;
        bl(ii_b).bbpim = [im_low_lims, im_low_lims + blob_p_sizes(ii_b) - 1];

        o.fprintf(repmat('\b', [1, num_chars]));
    end
    o.fprintf('Done in %g s\n', toc(t));
end