function shape_funcs = gtDefShapeFunctionsFwdProj(sampler, varargin)

    conf = struct( ...
        'shape_function_type', 'uvw', ...
        'recenter_sf', true, ...
        'data_type', 'double', ...
        'factor', 5, ...
        'rspace_oversampling', 1, ...
        'rspace_voxel_size', [1 1 1], ...
        'interpolated_voxel', false, ...
        'projector', 'nearest');
    conf = parse_pv_pairs(conf, varargin);

    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', conf.recenter_sf);
    fedpars = struct( ...
        'bltype', conf.data_type, ...
        'dcomps', [1 1 1 0 0 0 0 0 0] == 1, ...
        'defmethod', 'rod_rightstretch', ...
        'detector', repmat( fed_pars_detector, [num_detectors, 1]), ...
        'projector', conf.projector );

    voxel_size = sampler.stats.sampling.gaps;

    space_res = tand(sampler.estimate_maximum_resolution() ./ 2);
    voxel_sampling = max(ceil(voxel_size ./ space_res * conf.factor), 5);

    num_sub_voxels = prod(voxel_sampling);
    ones_nsv = ones(num_sub_voxels, 1);

    gv = struct('pcs', [], 'ind', [], 'd', [], 'pow', []);

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

    ref_gr = sampler.get_reference_grain();
    gv.pcs = compute_centers(ref_gr.center, conf.rspace_voxel_size, conf.rspace_oversampling)';

    % computing the subsampling, used to deterine the shape functions
    voxel_size = sampler.stats.sampling.gaps;
    if (conf.interpolated_voxel)
        half_voxel_size = voxel_size; % Enlarging to neighbours
        steps = voxel_size ./ voxel_sampling * 2;
    else
        half_voxel_size = voxel_size ./ 2;
        steps = voxel_size ./ voxel_sampling;
    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 (conf.interpolated_voxel)
        gv.pow = prod(1 - abs(gv.d ./ half_voxel_size(ones_nsv, :)'), 1);
    else
        gv.pow = ones(1, num_sub_voxels, conf.data_type);
    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 %s shape functions: ', upper(conf.shape_function_type))
    c = tic();
    for ii_g = 1:num_ors
        num_chars = fprintf('%03d/%03d', ii_g, num_ors);

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

        or_gv = gv;
        or_gv.d = gtMathsRodSum(gv.d, or.R_vector');

        switch (lower(conf.shape_function_type))
            case 'w'
                shape_funcs{ii_g} = gtDefFwdProjGvdm2W(or, sel_bls, or_gv, ...
                    fedpars, sampler.parameters, det_ind, false);
            case 'nw'
                eta_step = 2 * atand(space_res(1));
                shape_funcs{ii_g} = gtDefFwdProjGvdm2NW(or, sel_bls, or_gv, ...
                    fedpars, sampler.parameters, eta_step, det_ind, false);
            case 'uvw'
                shape_funcs{ii_g} = gtDefFwdProjGvdm(or, sel_bls, or_gv, ...
                    fedpars, sampler.parameters, det_ind, false);
        end

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

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

function centers = compute_centers(ref_center, rspace_voxel_size, rspace_oversampling)
    if (numel(rspace_oversampling) == 1)
        rspace_oversampling = rspace_oversampling([1 1 1]);
    end
    max_dists_from_center = rspace_voxel_size .* (1 - 1 ./ rspace_oversampling) / 2;
    tot_voxels = prod(rspace_oversampling);

    centers = zeros(tot_voxels, 3);

    counter = 1;

    for ss_z = linspace(-max_dists_from_center(3), max_dists_from_center(3), rspace_oversampling(3))
        for ss_y = linspace(-max_dists_from_center(2), max_dists_from_center(2), rspace_oversampling(2))
            for ss_x = linspace(-max_dists_from_center(1), max_dists_from_center(1), rspace_oversampling(1))

                centers(counter, :) = ref_center + [ss_x, ss_y, ss_z];

                counter = counter + 1;
            end
        end
    end
end