gtDefShapeFunctionsCreateW.m

function sf = gtDefShapeFunctionsCreateW(sampler)
% FUNCTION sf = gtDefShapeFunctionsCreateW(sampler)
%

    p = sampler.parameters;

    det_ind = sampler.detector_index;

    omstep = gtAcqGetOmegaStep(p, det_ind);
    omstep_rod = tand(omstep / 2);

    ref_gr = sampler.get_reference_grain();

    ref_ond = ref_gr.proj(det_ind).ondet;
    ref_inc = ref_gr.proj(det_ind).included;
    ref_sel = ref_gr.proj(det_ind).selected;

    ref_ab_sel = ref_ond(ref_inc(ref_sel));
    ref_ab_inc = ref_ond(ref_inc);

    num_inc_blobs = numel(ref_ab_sel);

    if (isfield(ref_gr.allblobs, 'plcry'))
        pl_cry = ref_gr.allblobs(det_ind).plcry(ref_ab_sel, :);
    else
        pl_samd = ref_gr.allblobs(det_ind).plorig(ref_ab_sel, :);
        g = gtMathsRod2OriMat(ref_gr.R_vector');
        pl_cry = gtVectorLab2Cryst(pl_samd, g)';
    end

    sub_space_size = tand(sampler.sampling_res / 2);
    ospace_samp_size = sampler.get_orientation_sampling_size();
    ospace_samp_size = ospace_samp_size{1};

    num_ospace_oversampling = prod(sampler.ss_factor);
    if (num_ospace_oversampling == 1)
        ospace_samp_lattice = sampler.lattice(det_ind).gr;
    else
        ospace_samp_lattice = cat(4, sampler.lattice_ss(det_ind).gr{:});
        ospace_samp_lattice = reshape(ospace_samp_lattice, [sampler.ss_factor, ospace_samp_size]);
        ospace_samp_lattice = permute(ospace_samp_lattice, [1 4 2 5 3 6]);

        ospace_samp_size = ospace_samp_size .* sampler.ss_factor;
        sub_space_size = sub_space_size ./ sampler.ss_factor;

        ospace_samp_lattice = reshape(ospace_samp_lattice, ospace_samp_size);
    end

    % Let's first compute the W extent
    half_rspace_sizes = sub_space_size / 2;

    fprintf('Computing bounds of W shape functions: ')
    c = tic();
    tot_orientations = numel(ospace_samp_lattice);

    oris_bounds_size = ospace_samp_size + 1;
    oris_lims_min = ospace_samp_lattice{1, 1, 1}.R_vector - half_rspace_sizes;
    oris_lims_max = ospace_samp_lattice{end, end, end}.R_vector + half_rspace_sizes;

    oris_steps_x = linspace(oris_lims_min(1), oris_lims_max(1), oris_bounds_size(1));
    oris_steps_y = linspace(oris_lims_min(2), oris_lims_max(2), oris_bounds_size(2));
    oris_steps_z = linspace(oris_lims_min(3), oris_lims_max(3), oris_bounds_size(3));

    oris_bounds = cell(oris_bounds_size);
    for dx = 1:oris_bounds_size(1)
        for dy = 1:oris_bounds_size(2)
            for dz = 1:oris_bounds_size(3)
                r_vec = [oris_steps_x(dx), oris_steps_y(dy), oris_steps_z(dz)];
                oris_bounds{dx, dy, dz} = struct( ...
                    'phaseid', ref_gr.phaseid, ...
                    'center', ref_gr.center, 'R_vector', r_vec );
            end
        end
    end

    chunk_size = 250;
    for ii = 1:chunk_size:numel(oris_bounds)
        end_chunk = min(ii+chunk_size, numel(oris_bounds));
        num_chars_gr_ii = fprintf('%05d/%05d', ii, numel(oris_bounds));

        oris_bounds(ii:end_chunk) = gtCalculateGrain_p( ...
            oris_bounds(ii:end_chunk), p, ...
            'ref_omind', ref_gr.allblobs(det_ind).omind, ...
            'included', ref_ab_inc );

        fprintf(repmat('\b', [1 num_chars_gr_ii]));
    end
    fprintf('Done in %g seconds.\n', toc(c))

    ref_ws = ref_gr.allblobs(det_ind).omega(ref_ab_sel) ./ omstep;

    fprintf('Precomputing shared values..')
    % Forming the grid in orientation space, on the plane that is
    % perpendicular to the axis of rotation (z-axis in the w case),
    % where to sample the sub-regions of the orientation-space voxels
    space_res = tand(sampler.estimate_maximum_resolution() ./ 2);
    num_poins = max(ceil(sub_space_size ./ space_res ./ 2), 5);
    pos_x = linspace(-half_rspace_sizes(1), half_rspace_sizes(1), num_poins(1)+1);
    pos_y = linspace(-half_rspace_sizes(2), half_rspace_sizes(2), num_poins(2)+1);

    renorm_coeff = (pos_x(2) - pos_x(1)) * (pos_y(2) - pos_y(1)) / prod(sub_space_size);

    pos_x = (pos_x(1:end-1) + pos_x(2:end)) / 2;
    pos_y = (pos_y(1:end-1) + pos_y(2:end)) / 2;

    pos_x_exp = reshape(pos_x, 1, [], 1);
    pos_x_exp = pos_x_exp(1, :, ones(1, num_poins(2)));
    pos_y_exp = reshape(pos_y, 1, 1, []);
    pos_y_exp = pos_y_exp(1, ones(1, num_poins(1)), :);
    pos_z_exp = zeros(1, num_poins(1), num_poins(2));

    tot_sampling_points = num_poins(1) * num_poins(2);
    ones_tot_sp = ones(tot_sampling_points, 1);

    % There is no point in computing the same rotation matrices each
    % time, so we pre compute them here and store them for later
    c = tic();
    % selecting Ws to find the interval
    ori_bounds = [oris_bounds{:}];
    allblobs = [ori_bounds(:).allblobs];
    w_tabs_int = [allblobs(:).omega];
    w_tabs_int = w_tabs_int(ref_sel, :) ./ omstep;
    w_tabs_int = gtGrainAnglesTabularFix360deg(w_tabs_int, ref_ws, p);
    w_tabs_int = reshape(w_tabs_int, [], oris_bounds_size(1), oris_bounds_size(2), oris_bounds_size(3));

    rotcomp_z = gtMathsRotationMatrixComp([0, 0, 1]', 'col');

    rotcomp_w = gtMathsRotationMatrixComp(p.labgeo.rotdir', 'col');

    % Precomputing A matrix components for all the used omegas
    ref_round_int_ws = round(ref_ws + 0.5) - 0.5;
    ref_round_ws = ref_round_int_ws .* omstep;
    rot_w = gtMathsRotationTensor(ref_round_ws, rotcomp_w);

    beamdirs = p.labgeo.beamdir * reshape(rot_w, 3, []);
    beamdirs = reshape(beamdirs, 3, [])';

    Ac_part = beamdirs * rotcomp_z.cos;
    As_part = beamdirs * rotcomp_z.sin;
    Acc_part = beamdirs * rotcomp_z.const;

    Ac_part = Ac_part(:, :, ones_tot_sp);
    As_part = As_part(:, :, ones_tot_sp);
    Acc_part = Acc_part(:, :, ones_tot_sp);

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

    sf = cell(tot_orientations, 1);

    [o_x, o_y, o_z] = ind2sub(oris_bounds_size-1, 1:tot_orientations);

    fprintf('Computing W shape functions: ')
    c = tic();
    for ii_g = 1:tot_orientations
        if (mod(ii_g, chunk_size) == 1)
            num_chars = fprintf('%05d/%05d', ii_g, tot_orientations);
        end

        or = ospace_samp_lattice{ii_g};

        % Extracting ospace boundaries for the given voxel
        w_tab_int = w_tabs_int(:, o_x(ii_g):o_x(ii_g)+1, o_y(ii_g):o_y(ii_g)+1, o_z(ii_g):o_z(ii_g)+1);
        w_tab_int = reshape(w_tab_int, [], 8);

        lims_w_proj_g_int = round([min(w_tab_int, [], 2), max(w_tab_int, [], 2)]);

        % Forming the grid in orientation space, on the plane that is
        % perpendicular to the axis of rotation (z-axis in the w case),
        % where to sample the sub-regions of the orientation-space voxel
        % that we want to integrate.
        r_vecs = [pos_x_exp; pos_y_exp; pos_z_exp];
        r_vecs = reshape(r_vecs, 3, []);
        r_vecs = or.R_vector(ones_tot_sp, :)' + r_vecs;

        gs = gtMathsRod2OriMat(r_vecs);
        % Unrolling and transposing the matrices at the same time
        % because we need g^-1
        gs = reshape(gs, 3, [])';

        ys = gs * pl_cry;
        ys = reshape(ys, 3, [], num_inc_blobs);
        ys = permute(ys, [3 1 2]);

        % Starting here to compute initial shifts aligned with the
        % rounded omega, by frst doing three scalar products
        Ac = sum(Ac_part .* ys, 2);
        As = sum(As_part .* ys, 2);
        Acc = sum(Acc_part .* ys, 2);

        Ac = reshape(Ac, num_inc_blobs, tot_sampling_points);
        As = reshape(As, num_inc_blobs, tot_sampling_points);
        Acc = reshape(Acc, num_inc_blobs, tot_sampling_points);

        D = Ac .^ 2 + As .^ 2;

        % Precomputing useful values, like sinthetas
        ominds = or.allblobs(det_ind).omind(ref_sel);
        ssp = ((ominds == 1) | (ominds == 2));
        ss  = ssp - ~ssp;
        sinths = ss .* or.allblobs(det_ind).sintheta(ref_sel);
        sinths = sinths(:, ones(tot_sampling_points, 1));

        CC = Acc + sinths;
        DD = D - CC .^ 2;
        E  = sqrt(DD);
        ok = DD > 0;

        % Apparently it is inverted, compared to the omega prediction
        % function. Should be investigated
        ssp = ~((ominds == 1) | (ominds == 3));
        ss  = ssp - ~ssp;
        ss = ss(:, ones(tot_sampling_points, 1));

        % Shift along z in orientation space, to align with the images
        % in proections space
        dz = (-As + ss .* ok .* E) ./ (CC - Ac);

        ws_disps = [ ...
            (lims_w_proj_g_int(:, 1) - 0.5 - ref_round_int_ws), ...
            (lims_w_proj_g_int(:, 2) + 0.5 - ref_round_int_ws) ];
        ws_disps = tand(ws_disps .* omstep / 2);

        num_ds = lims_w_proj_g_int(:, 2) - lims_w_proj_g_int(:, 1) + 2;

        bbwims = cell(num_inc_blobs, 1);
        intms = cell(num_inc_blobs, 1);
        bbsizes = cell(num_inc_blobs, 1);

        % This is the most expensive operation because the different
        % sizes of the omega spread, which depends on the blob itself,
        % don't allow to have vectorized operations on all the blobs at
        % the same time
        for ii_b = 1:num_inc_blobs
            d = (ws_disps(ii_b, 1):omstep_rod:ws_disps(ii_b, 2))';
            d = d(:, ones(tot_sampling_points, 1)) + dz(ii_b * ones(num_ds(ii_b), 1), :);

            upper_limits = min(d(2:end, :), half_rspace_sizes(3));
            lower_limits = max(d(1:end-1, :), -half_rspace_sizes(3));

            % integration is simple and linear: for every w and
            % xy-position, we consider the corresponding segment on the
            % z-direction
            ints = upper_limits - lower_limits;
            ints(ints < 0) = 0;
            ints = sum(ints, 2);

            bbwims{ii_b} = lims_w_proj_g_int(ii_b, :);
            intms{ii_b} = renorm_coeff * ints;
            bbsizes{ii_b} = numel(ints);
        end

        sf_bls_w = gtFwdSimBlobDefinition('sf_w', num_inc_blobs);
        [sf_bls_w(:).bbwim] = deal(bbwims{:});
        [sf_bls_w(:).intm] = deal(intms{:});
        [sf_bls_w(:).bbsize] = deal(bbsizes{:});
        sf{ii_g} = sf_bls_w;

        if (mod(ii_g, chunk_size) == 1)
            fprintf(repmat('\b', [1, num_chars]));
        end
    end
    fprintf('Done in %g seconds.\n', toc(c))

    if (num_ospace_oversampling > 1)
        fprintf('Re-grouping oversampled W shape functions: ')
        c = tic();
        ospace_samp_size = sampler.get_orientation_sampling_size();
        ospace_samp_size = ospace_samp_size{1};
        tot_orient = sampler.get_total_num_orientations();

        regroup = [ ...
            sampler.ss_factor(1), ospace_samp_size(1), ...
            sampler.ss_factor(2), ospace_samp_size(2), ...
            sampler.ss_factor(3), ospace_samp_size(3) ];
        sf = reshape(sf, regroup);
        sf = permute(sf, [1 3 5 2 4 6]);
        sf = reshape(sf, num_ospace_oversampling, tot_orient);
        fprintf('Done in %g seconds.\n', toc(c))
    end
end