Added some more bits from Lorenzo's work: the active volume generation

Signed-off-by: Nicola Vigano <nicola.vigano@esrf.fr>

5_reconstruction/gtAstraBackproject.m

+function vol = gtAstraBackproject(proj_stack, geometry, vol_size)
+
+    usize = size(proj_stack, 1);
+    vsize = size(proj_stack, 3);
+
+    % projection geometry (type, #rows, #columns, vectors)
+    proj_geom = astra_create_proj_geom('parallel3d_vec', vsize, usize, geometry);
+
+    % volume geometry (x, y, z)
+    vol_geom = astra_create_vol_geom(vol_size(2), vol_size(1), vol_size(3));
+
+    astra_projector_id = astra_create_projector('cuda3d', proj_geom, vol_geom);
+
+    vol = astra_mex_direct('BP3D', astra_projector_id, proj_stack);
+
+    astra_mex_projector('delete', astra_projector_id);
+end

7_fed2/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));
+
+    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, :) = true;
+    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
+