gtDefShapeFunctionsFwdProjUVW.m

function shape_funcs = gtDefShapeFunctionsFwdProjUVW(sampler, factor, dtype, recenter_sf)

    if (~exist('dtype', 'var') || isempty(dtype))
        dtype = 'single';
    end

    if (~exist('recenter_sf', 'var') || isempty(recenter_sf))
        recenter_sf = true;
    end

    interpolated_voxel = false;

    num_detectors = numel(sampler.parameters.detgeo);
    det_ind = sampler.detector_index;

    labgeo = sampler.parameters.labgeo;
    labgeo.rotcomp = gtMathsRotationMatrixComp(labgeo.rotdir', 'col');
    sampler.parameters.labgeo = labgeo;

    fed_pars_detector = struct(...
        'blobsizeadd', [0 0 0], 'psf', [], 'apply_uv_shift', recenter_sf);
    fedpars = struct( ...
        'bltype', dtype, ...
        'dcomps', [1 1 1 0 0 0 0 0 0] == 1, ...
        'defmethod', 'rod_rightstretch', ...
        'detector', repmat( fed_pars_detector, [num_detectors, 1]));

    voxel_factor = factor([1 1 1]);
    num_sub_voxels = prod(voxel_factor);
    ones_nsv = ones(num_sub_voxels, 1);

    gv = struct('cs', [], 'ind', [], 'dm', [], 'pow', []);

    [gv.ind(:, 1), gv.ind(:, 2), gv.ind(:, 3)] = ind2sub(voxel_factor, 1:num_sub_voxels);
    gv.ind = gv.ind';

    ref_gr = sampler.get_reference_grain();
    gv.cs = ref_gr.center(ones_nsv, :)';

    % computing the subsampling, used to deterine the shape functions
    voxel_size = sampler.stats.sampling.gaps;
    if (interpolated_voxel)
        half_voxel_size = voxel_size; % Enlarging to neighbours
        steps = voxel_size ./ voxel_factor * 2;
    else
        half_voxel_size = voxel_size ./ 2;
        steps = voxel_size ./ voxel_factor;
    end
    half_steps = steps ./ 2;
    beg_pos = - half_voxel_size + half_steps;
    end_pos = + half_voxel_size - half_steps;

    x_steps = beg_pos(1):steps(1):end_pos(1);
    y_steps = beg_pos(2):steps(2):end_pos(2);
    z_steps = beg_pos(3):steps(3):end_pos(3);

    num_steps = [numel(x_steps), numel(y_steps), numel(z_steps)];
    ones_x_steps = ones(1, num_steps(1));

    gv.d = zeros(3, prod(num_steps));

    pos = 0;
    for ii_z = z_steps
        for ii_y = y_steps
            gv.d(:, (pos+1):(pos+num_steps(1))) ...
                = [ x_steps; ones_x_steps * ii_y; ones_x_steps * ii_z];

            pos = pos + num_steps(1);
        end
    end

    if (interpolated_voxel)
        gv.pow = prod(1 - abs(gv.d ./ half_voxel_size(ones_nsv, :)'), 1);
    else
        gv.pow = ones(1, num_sub_voxels, dtype);
    end
    gv.pow = gv.pow / sum(gv.pow);

    gv.used_ind = find(gv.pow > 0);

    sel_bls = find(sampler.selected);

    % Only one lattice allowed for the moment!
    num_ors = numel(sampler.lattice(1).gr);
    shape_funcs = cell(num_ors, 1);

    fprintf('Computing UVW shape functions: ')
    c = tic();
    for ii_g = 1:num_ors
        num_chars = fprintf('%03d/%03d', ii_g, num_ors);

        or = sampler.lattice(1).gr{ii_g};

        shape_funcs{ii_g} = gtDefFwdProjGvdm(or, sel_bls, gv, ...
            fedpars, sampler.parameters, det_ind, false);

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

    fprintf('Done in %g seconds.\n', toc(c))
end