gtDefFwdProjGvdm.m

function bl = gtDefFwdProjGvdm(grain, ref_sel, 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)

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

    detgeo = parameters.detgeo(det_ind);
    labgeo = parameters.labgeo;
    samgeo = parameters.samgeo;

    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);
    or_uvw = grain.allblobs.detector(det_ind).uvw(ref_sel, :);
    ref_ws = or_uvw(:, 3);

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

    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)
        uvw_shifts = [or_uvw(:, 1:2), round(or_uvw(:, 3))]';
    else
        uvw_shifts = round(or_uvw)';
    end

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

    num_oversampling = size(gv.d, 3);
    gvpow = 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', 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

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

        uvw = cell(nbl, 1);
        uvw_min = zeros(nbl, 3);
        uvw_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)
            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, pllab, ~, rot_l2s] = 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

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

            rot_s2l = permute(rot_l2s, [2 1 3]);
            gvcs_lab = gtGeoSam2Lab(gvpcs', rot_s2l, labgeo, samgeo, false, [], element_wise_gcs);

            uvw_bl = gtFedPredictUVWMultiple([], dvec_lab, gvcs_lab', ...
                detgeo.detrefpos', detgeo.detnorm', detgeo.Qdet, ...
                [detgeo.detrefu, detgeo.detrefv]', om, labgeo.omstep);

            uvw_bl(3, :) = gtGrainAnglesTabularFix360deg(uvw_bl(3, :), ...
                ref_ws(ii_b), parameters);

            % Transforming into offsets from 
            uvw_bl = uvw_bl - uvw_shifts(:, ii_b * ones(1, nv));

            uvw{ii_b} = uvw_bl;

            uvw_min(ii_b, :) = min(uvw{ii_b}, [], 2)';
            uvw_max(ii_b, :) = max(uvw{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')
        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(uvw_min(:, 1:2), [], 1)), abs(max(uvw_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_w_sizes = ceil(uvw_max(:, 3)) - floor(uvw_min(:, 3)) + blob_size_add(3) + 1;

        blob_orig_w_shifts = 1 - floor(uvw_min(:, 3));

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

        tmp_bl = gtFwdSimBlobDefinition('blob', 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_uvw_shift = [blob_orig_uv_shift, blob_orig_w_shifts(ii_b)]';
            blob_uvw_size = [blob_uv_size blob_w_sizes(ii_b)];

            % Detector coordinates U,V in blob
            uvw_bl = uvw{ii_b} + blob_orig_uvw_shift(:, ones(1, nv));

            % Let's now filter valid voxels
            uvw_bl = uvw_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(uvw_bl, ...
                    gvpow_v, fedpars.bltype);
            else
                ucbl8 = round(uvw_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_uvw_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_uvw_size)
                disp('Min projected U,V,W coordinates:')
                fprintf('\t%g\t%g\t%g \t(rounded: %d %d %d )\n', ...
                    min(uvw_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(uvw_bl, [], 1), max(ucbl8, [], 1))
                new_exc = GtFedExceptionFwdProj(mexc, det_ind, [0 0 0]);
                throw(new_exc)
            end

            tmp_bl(ii_b).bbsize = blob_uvw_size;

            im_low_lims = uvw_shifts(:, ii_b) - blob_orig_uvw_shift + 1;
            tmp_bl(ii_b).bbuim = [im_low_lims(1), im_low_lims(1) + blob_uvw_size(1) - 1];
            tmp_bl(ii_b).bbvim = [im_low_lims(2), im_low_lims(2) + blob_uvw_size(2) - 1];
            tmp_bl(ii_b).bbwim = [im_low_lims(3), im_low_lims(3) + blob_uvw_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(:).bbwim);
            tmp_bl_bbs = reshape(tmp_bl_bbs, nbl, 3, 2);

            bl_bbs = cat(1, bl(:).bbuim, bl(:).bbvim, bl(:).bbwim);
            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');

                    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.bbwim = 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;

        wproj = sum(sum(abs(bl(ii_b).intm), 1), 2);
        bl(ii_b).mbbw = [find(wproj, 1, 'first'), find(wproj, 1, 'last')] + bl(ii_b).bbwim(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).mbbw(2) - bl(ii_b).mbbw(1) + 1 ];
    end
    o.fprintf('\b\b: Done in %g s\n', toc(t));
end