gtFedActiveVolumeVoxels.m

function active_vol = gtFedActiveVolumeVoxels(phase_vol, parameters_inline, varargin)
% GTFEDACTIVEVOLUMEVOXELS    Backprojects the absorption
% image from a mask on the sample volume, obtaining the diffracting
% voxels for a given diffraction geometry.
%
% NOTE: This function was written for the Mg monocrystal. If you need to
% change it, for another dataset, please make it ust more customizable!
%
% gtFedActiveVolumeVoxels(phase_vol)
% ---------------------------------------------------------------------
%
% INPUT
%   phase_vol         : assembled phase volume reconstructed with 6D
%                       algorithm
%   parameters_inline : parameters for an in-line acquisition
%
% OUTPUT
%   active_vol        : three volumes with modulated intensities
%                       according to scattering powers, one for each
%                       considered reflection
%
% OPTIONAL INPUT
%   mask              : {'pencil'} | 'line'
%   det_index         : {1}
%   vert_slice        : {[]}, and empty field selects the central slice
%   pencil_spacing    : {0.05}, determines the spacing between each window
%                       and it is in microns

    conf = struct(...
        'mask', 'pencil', ...
        'det_index', 1, ...
        'vert_slice', [], ...
        'pencil_spacing', 0.050 );
    conf = parse_pv_pairs(conf, varargin);

    if (~exist('parameters_inline', 'var') || isempty(parameters_inline))
        parameters_inline = gtLoadParameters();
    end

    detgeo = parameters_inline.detgeo(conf.det_index);
    labgeo = parameters_inline.labgeo;
    samgeo = parameters_inline.samgeo;
    recgeo = parameters_inline.recgeo(conf.det_index);
    acq = parameters_inline.acq(conf.det_index);

    vol_size_xyz = size(phase_vol.intvol);

    det_mask = zeros(vol_size_xyz(2), vol_size_xyz(3), 'single');

    if (isempty(conf.vert_slice))
        conf.vert_slice = round(vol_size_xyz(2) / 2);
    end

    fprintf('Creating a %s mask..', conf.mask)
    c = tic();
    switch(conf.mask)
        case 'pencil'
            % Creation of the pencil mask: 7 um windows with 50 um period,
            % 1 vx vertical thickness
            vox_size = recgeo.voxsize(3);
            coords = ( vox_size:conf.pencil_spacing:(vol_size_xyz(3) * vox_size) ) / vox_size;
            [inds, ints] = gtMathsGetInterpolationIndices(coords', ones(numel(coords), 1));

            good_inds = ints > eps('single');
            inds = inds(good_inds);
            ints = ints(good_inds);
            v = accumarray(inds, ints, [vol_size_xyz(3), 1])';

            det_mask(conf.vert_slice, :) = v;
        case 'line'
            % Creation of the line mask: 1 vx vertical thickness
            det_mask(conf.vert_slice, :) = 1;
    end
    fprintf('\b\b: Done in %g seconds.\n', toc(c))

    % Manual correction of intensity volume
    phase_vol_corr = single(phase_vol.intvol);

    % Back-projection of absorption image into sample volume
    pmos = [-88.2 91.8 -26.2];
    omegas = pmos + 90;
    num_rots = numel(pmos);

    proj_geoms = gtGeoProjForReconstruction([], omegas, [], acq.bb, ...
        zeros(num_rots, 2), detgeo, labgeo, samgeo, recgeo, ...
        'ASTRA_absorption');

    active_vol = zeros([vol_size_xyz, num_rots], 'single');

    for irefl = 1:num_rots
        fprintf('Backprojecting, angle = %g, pmo = %g..', ...
            omegas(irefl), pmos(irefl))
        c = tic();

        proj_stack(:, 1, :) = det_mask;
        vol = gtAstraBackproject(proj_stack, proj_geoms(irefl, :), vol_size_xyz);

        active_vol(:, :, :, irefl) = vol .* phase_vol_corr;
        fprintf('\b\b: Done in %g seconds.\n', toc(c))
    end
end