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;

    det_ind = sampler.detector_index;

    labgeo = sampler.parameters.labgeo;
%     samgeo = sampler.parameters.samgeo;
%     recgeo = sampler.parameters.recgeo(det_ind);
%     detgeo = sampler.parameters.detgeo(det_ind);

    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( ...
        'dcomps', [1 1 1 0 0 0 0 0 0] == 1, ...
        'defmethod', 'rod_rightstretch', ...
        'detector', repmat( fed_pars_detector, [numel(sampler.parameters.detgeo), 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);
    nuv = numel(gv.used_ind);
    ones_nuv = ones(1, nuv);

    sel_bls = find(sampler.selected);
    inc_bls = sampler.ondet(sampler.included(sampler.selected));

    % computing the size of the blobs
    shape_funcs_size(:, 3) = ceil(2 * atand(sampler.stats.sampling.gaps(3)) ...
        / labgeo.omstep .* ref_gr.allblobs.lorentzfac(inc_bls)) + 1;
    shape_funcs_size(:, 1) = sampler.stats.proj.max_pixel_dist_u;
    shape_funcs_size(:, 2) = sampler.stats.proj.max_pixel_dist_v;
    if (interpolated_voxel)
        shape_funcs_size = shape_funcs_size * 2;
    end
    shape_funcs_size = 2 * ceil(shape_funcs_size) + 1; % Imposing odd shape! -> for centering
    half_sf_size = floor(shape_funcs_size / 2);

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

    nbl = numel(sel_bls);

    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};

        bl = repmat(gtFedBlobDefinition(), nbl, 1);

        or_uvw = or.allblobs.detector(det_ind).uvw(sel_bls, :);

        for ii_b = 1:nbl
            % Load diffraction paramaters of blobs

            % Initial strained non-signed plane normal in SAMPLE reference
            bl(ii_b).plorig = or.allblobs.plorig(sel_bls(ii_b), :);
            bl(ii_b).sinth0 = or.allblobs.sintheta(sel_bls(ii_b));  % initial sin(theta)
            bl(ii_b).omega0 = or.allblobs.omega(sel_bls(ii_b));     % initial omega
            bl(ii_b).omind  = or.allblobs.omind(sel_bls(ii_b));     % initial omega index (1..4)

            % Initial centroid in image
            bl(ii_b).u0im = or_uvw(ii_b, :);

            % Projection information!!
            bl(ii_b).bbuim = round(or_uvw(ii_b, 1) + [-half_sf_size(ii_b, 1), half_sf_size(ii_b, 1)]);
            bl(ii_b).bbvim = round(or_uvw(ii_b, 2) + [-half_sf_size(ii_b, 2), half_sf_size(ii_b, 2)]);
            bl(ii_b).bbwim = round(or_uvw(ii_b, 3) + [-half_sf_size(ii_b, 3), half_sf_size(ii_b, 3)]);
            bl(ii_b).bbsize = [ ...
                (bl(ii_b).bbuim(2) - bl(ii_b).bbuim(1) + 1), ...
                (bl(ii_b).bbvim(2) - bl(ii_b).bbvim(1) + 1), ...
                (bl(ii_b).bbwim(2) - bl(ii_b).bbwim(1) + 1) ];

            % Blob origin and center [u, v, w]:
            bl(ii_b).bbor   = [bl(ii_b).bbuim(1)-1, bl(ii_b).bbvim(1)-1, bl(ii_b).bbwim(1)-1];

            bl(ii_b).bbcent = [ ...
                round( (bl(ii_b).bbuim(1) + bl(ii_b).bbuim(2) ) / 2 ), ...
                round( (bl(ii_b).bbvim(1) + bl(ii_b).bbvim(2) ) / 2 ), ...
                round( (bl(ii_b).bbwim(1) + bl(ii_b).bbwim(2) ) / 2 ) ];

            % u in blob = u in image minus u at blob origin
            bl(ii_b).u0bl = bl(ii_b).u0im - bl(ii_b).bbor;
            bl(ii_b).uv_shift = bl(ii_b).u0im(1:2) - bl(ii_b).bbcent(1:2);

            bl(ii_b).int = zeros(bl(ii_b).bbsize, dtype);
            bl(ii_b).intm = [];
        end

        for ii_b = nbl:-1:1
            gvb(ii_b).uc      = zeros(3, nuv, dtype);
            gvb(ii_b).ucd     = zeros(3, nuv, dtype);
            gvb(ii_b).ucblmod = zeros(3, nuv, dtype);
            gvb(ii_b).ucbl8   = zeros(1, nuv*8, dtype);
            gvb(ii_b).ucbl8ini= zeros(1, nuv*8, dtype);
            gvb(ii_b).ucbl8ok = 0;

            % Initial u of gv centre relative to blob origin
            gvb(ii_b).uc0bl = cast(bl(ii_b).u0im - bl(ii_b).bbor, dtype)';
            gvb(ii_b).uc0bl = gvb(ii_b).uc0bl(:, ones_nuv);
        end

        [gvb, bl] = gtFedFwdProjExact(gv, gvb, bl, fedpars, sampler.parameters, det_ind, false);

        % Cutting borders
        for ii_b = 1:nbl
            diff_inds_w = [...
                find(sum(sum(bl(ii_b).int, 1), 2), 1, 'first') - 1, ...
                find(sum(sum(bl(ii_b).int, 1), 2), 1, 'last') - bl(ii_b).bbsize(3); ...
                ];
            bl(ii_b).bbwim = bl(ii_b).bbwim + diff_inds_w;

            bl(ii_b).mbbu = bl(ii_b).bbuim;
            bl(ii_b).mbbv = bl(ii_b).bbvim;
            bl(ii_b).mbbw = bl(ii_b).bbwim;

            inds_w = [1, bl(ii_b).bbsize(3)] + diff_inds_w;

            bl(ii_b).intm = bl(ii_b).int(:, :, inds_w(1):inds_w(2));
            bl(ii_b).int = [];

            bl(ii_b).bbsize(3) = size(bl(ii_b).intm, 3);
        end

        shape_funcs{ii_g} = bl;

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

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