function bl = gtDefFwdProjGvdm2UVP(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)

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

    uinds = gv.used_ind;
    nv = numel(uinds);
    ones_bl = ones(nbl, 1);

    detgeo = parameters.detgeo(det_ind);
    labgeo = parameters.labgeo;
    samgeo = parameters.samgeo;
    diff_acq = parameters.diffractometer(det_ind);
    diff_ref = parameters.diffractometer(1);

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

    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

    % 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

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

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

    if (isfield(fedpars, 'detector') ...
            && isfield(fedpars.detector, 'psf') ...
            && ~isempty(fedpars.detector(det_ind).psf))
        psf = fedpars.detector(det_ind).psf;
        if (~iscell(psf))
            psf = { psf };
        end
        if (numel(psf) == 1)
            psf = psf(ones_bl);
        end
    else
        psf = {};
    end

    if (isfield(fedpars, 'detector') ...
            && isfield(fedpars.detector, 'apply_uv_shift') ...
            && ~isempty(fedpars.detector(det_ind).apply_uv_shift) ...
            && fedpars.detector(det_ind).apply_uv_shift)
        uvpw_shifts = [or_uvpw(:, 1:2), round(or_uvpw(:, 3))]';
    else
        uvpw_shifts = round(or_uvpw)';
    end

    linear_interp = ~(isfield(fedpars, 'projector') ...
        && isstruct(fedpars.projector) && isfield(fedpars.projector, 'interp') ...
        && strcmpi(fedpars.projector.interp, 'nearest'));

    num_oversampling = size(gv.d, 3);
    gvpow_uinds = gv.pow(1, uinds, :);

    % This usually happens for shape functions, where we only have one
    % center
    element_wise_gcs = size(gv.pcs, 2) ~= 1;

    bl = gtFwdSimBlobDefinition('blob_topo', nbl);

    for ii_ss = 1:num_oversampling
        if (element_wise_gcs)
            gvpcs = gv.pcs(:, uinds, ii_ss)';
        else
            gvpcs = gv.pcs(:, 1, ii_ss)';
        end
        if (size(gv.pow, 3) > 1)
            gvpow = gvpow_uinds(1, :, ii_ss);
        else
            gvpow = gvpow_uinds;
        end

        gvd = gv.d(:, uinds, ii_ss);

        % Deformation tensor (relative to reference state) from its components
        defT = gtFedDefTensorFromComps(gvd, fedpars.dcomps, fedpars.defmethod, 0);

        %%% Computation of indices
        o.fprintf(' - Super sampling %03d/%03d:\n   * Computing indices and bbsizes: ', ...
            ii_ss, num_oversampling)
        t = tic();

        uvp = cell(nbl, 1);
        uvp_min = zeros(nbl, 3);
        uvp_max = zeros(nbl, 3);

        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, pllab] = 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

            % Diffraction vector
            dvec_lab = gtFedPredictDiffVecMultiple(pllab, labgeo.beamdir');

            rot_s2l_tot = gtGeoDiffractometerTensor(diff_acq, 'sam2lab', ...
                'reference_diffractometer', diff_ref, ...
                'angles_rotation', omega(ii_b), ...
                'angles_basetilt', ph);

            gvcs_lab = gtGeoSam2Lab(gvpcs, rot_s2l_tot, labgeo, samgeo, ...
                false, 'element_wise', element_wise_gcs);

            uv_bl = gtFedPredictUVMultiple([], dvec_lab, gvcs_lab', ...
                detgeo.detrefpos', detgeo.detnorm', detgeo.Qdet, ...
                [detgeo.detrefu, detgeo.detrefv]');
            uvp_bl = [uv_bl; ph ./ phstep];

            % Transforming into offsets from 
            uvp_bl = bsxfun(@minus, uvp_bl, uvpw_shifts(:, ii_b));

            uvp{ii_b} = uvp_bl;

            uvp_min(ii_b, :) = min(uvp{ii_b}, [], 2)';
            uvp_max(ii_b, :) = max(uvp{ii_b}, [], 2)';

            if (any(uvp_min(ii_b, :) > 0) || any(uvp_max(ii_b, :) < 0))
                warning('gtDefFwdProjGvdm2UVP: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')
        psf_size = zeros(nbl, 2);
        if (~isempty(psf))
            for ii_b = 1:nbl
                psf_size(ii_b, :) = (size(psf{ii_b}) - 1) / 2;
            end
        end
        if (isfield(fedpars, 'detector'))
            blob_size_add = fedpars.detector(det_ind).blobsizeadd;
        else
            blob_size_add = fedpars.blobsizeadd(det_ind, :);
        end
        blob_size_add = blob_size_add + [max(psf_size, [], 1), 0];

        blob_uv_size = max(abs(min(uvp_min(:, 1:2), [], 1)), abs(max(uvp_max(:, 1:2), [], 1)));
        blob_uv_size = 2 * (ceil(blob_uv_size) + blob_size_add(1:2)) + 1;

        blob_orig_uv_shift = (blob_uv_size - 1) / 2 + 1;

        blob_p_sizes = ceil(uvp_max(:, 3)) - floor(uvp_min(:, 3)) + blob_size_add(3) + 1;

        blob_orig_p_shifts = 1 - floor(uvp_min(:, 3));

        o.fprintf('     Computed Blob UV size: [%d, %d]\n', blob_uv_size);
        o.fprintf('     Computed Blob Phi sizes: [%s]\n', sprintf(' %d', blob_p_sizes));
        o.fprintf('     Blob size add: [%d, %d, %d] (psf: [%d, %d])\n', ...
            blob_size_add, max(psf_size, [], 1));

        tmp_bl = gtFwdSimBlobDefinition('blob_topo', nbl);

        %%% Blob projection
        o.fprintf('   * Projecting volumes: ')
        t = tic();
        for ii_b = 1:nbl
            num_chars = o.fprintf('%03d/%03d', ii_b, nbl);

            blob_orig_uvp_shift = [blob_orig_uv_shift, blob_orig_p_shifts(ii_b)]';
            blob_uvp_size = [blob_uv_size blob_p_sizes(ii_b)];

            % Detector coordinates U,V in blob
            uvp_bl = uvp{ii_b} + blob_orig_uvp_shift(:, ones(1, nv));

            % Let's now filter valid voxels
            uvp_bl = uvp_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(uvp_bl, ...
                    gvpow_v, fedpars.bltype);
            else
                ucbl8 = round(uvp_bl);
                ints8 = gvpow_v;
            end

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

            tmp_bl(ii_b).bbsize = blob_uvp_size;

            im_low_lims = uvpw_shifts(1:3, ii_b) - blob_orig_uvp_shift + 1;
            tmp_bl(ii_b).bbuim = [im_low_lims(1), im_low_lims(1) + blob_uvp_size(1) - 1];
            tmp_bl(ii_b).bbvim = [im_low_lims(2), im_low_lims(2) + blob_uvp_size(2) - 1];
            tmp_bl(ii_b).bbpim = [im_low_lims(3), im_low_lims(3) + blob_uvp_size(3) - 1];

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

        if (num_oversampling == 1 || ii_ss == 1)
            bl = tmp_bl;
        else
            o.fprintf('   * Merging super-sampling blobs: ')
            t = tic();
            tmp_bl_bbs = cat(1, tmp_bl(:).bbuim, tmp_bl(:).bbvim, tmp_bl(:).bbpim);
            tmp_bl_bbs = reshape(tmp_bl_bbs, nbl, 3, 2);

            bl_bbs = cat(1, bl(:).bbuim, bl(:).bbvim, bl(:).bbpim);
            bl_bbs = reshape(bl_bbs, nbl, 3, 2);

            should_reallocate = any(any(tmp_bl_bbs ~= bl_bbs, 2), 3);

            container_blobs_bb = cat(3, ...
                min(bl_bbs(:, :, 1), tmp_bl_bbs(:, :, 1)), ...
                max(bl_bbs(:, :, 2), tmp_bl_bbs(:, :, 2)) ...
                );
            container_blobs_size = container_blobs_bb(:, :, 2) - container_blobs_bb(:, :, 1) + 1;

            shift_bls = bl_bbs(:, :, 1) - container_blobs_bb(:, :, 1);
            shift_tmp_bls = tmp_bl_bbs(:, :, 1) - container_blobs_bb(:, :, 1);

            for ii_b = 1:nbl
                if (should_reallocate(ii_b))
                    new_bl = gtFwdSimBlobDefinition('blob_topo');

                    new_bl.intm = zeros(container_blobs_size(ii_b, :), fedpars.bltype);
                    new_bl.intm = gtPlaceSubVolume(new_bl.intm, bl(ii_b).intm, shift_bls(ii_b, :));
                    new_bl.intm = gtPlaceSubVolume(new_bl.intm, tmp_bl(ii_b).intm, shift_tmp_bls(ii_b, :), [], 'summed');

                    new_bl.bbsize = container_blobs_size(ii_b, :);

                    new_bl.bbuim = reshape(container_blobs_bb(ii_b, 1, :), 1, []);
                    new_bl.bbvim = reshape(container_blobs_bb(ii_b, 2, :), 1, []);
                    new_bl.bbpim = reshape(container_blobs_bb(ii_b, 3, :), 1, []);

                    bl(ii_b) = new_bl;
                else
                    bl(ii_b).intm = bl(ii_b).intm + tmp_bl(ii_b).intm;
                end
            end
            o.fprintf('Done in %g s\n', toc(t));
        end
    end

    if (num_oversampling > 1)
        o.fprintf(' - Renormalizing to initial intensity:');
        t = tic();

        for ii_b = 1:nbl
            bl(ii_b).intm = bl(ii_b).intm / num_oversampling;

            o.fprintf(repmat('\b', [1, num_chars]));
        end

        o.fprintf('Done in %g s\n', toc(t));
    end

    if (~isempty(psf))
        o.fprintf(' - Applying PSF:');
        t = tic();

        for ii_b = 1:nbl
            num_chars = o.fprintf('%03d/%03d', ii_b, nbl);

            bl(ii_b).intm = convn(bl(ii_b).intm, psf{ii_b}, 'same');

            o.fprintf(repmat('\b', [1, num_chars]));
        end

        o.fprintf('Done in %g s\n', toc(t));
    end

    o.fprintf(' - Computing boundaries of measured blobs inside blob bounding box..')
    t = tic();
    for ii_b = 1:nbl
        uproj = sum(sum(abs(bl(ii_b).intm), 2), 3);
        bl(ii_b).mbbu = [find(uproj, 1, 'first'), find(uproj, 1, 'last')] + bl(ii_b).bbuim(1) - 1;

        vproj = sum(sum(abs(bl(ii_b).intm), 1), 3);
        bl(ii_b).mbbv = [find(vproj, 1, 'first'), find(vproj, 1, 'last')] + bl(ii_b).bbvim(1) - 1;

        pproj = sum(sum(abs(bl(ii_b).intm), 1), 2);
        bl(ii_b).mbbp = [find(pproj, 1, 'first'), find(pproj, 1, 'last')] + bl(ii_b).bbpim(1) - 1;

        bl(ii_b).mbbsize = [bl(ii_b).mbbu(2) - bl(ii_b).mbbu(1) + 1, ...
                          bl(ii_b).mbbv(2) - bl(ii_b).mbbv(1) + 1, ...
                          bl(ii_b).mbbp(2) - bl(ii_b).mbbp(1) + 1 ];
    end
    o.fprintf('\b\b: Done in %g s\n', toc(t));

    for ii_b = 1:nbl
        bl(ii_b).intensity = sum(bl(ii_b).intm(:));
    end
end