Creation of pencil mask for pencil beam illumination

Signed-off-by: Lorenzo Valzania <valzania@rnice7-0308.esrf.fr>

7_fed2/gtFedCreateGrainPars.m

+function [fedpars, dmvol, intvol] = gtFedCreateGrainPars(grain, grain_rec, parameters, dmvol, intvol)
+% [fedpars, dmvol, intvol] = gtFedCreateGrainPars(grain, parameters, dmvol, intvol)
+    
+    %%%%%%%%%%%%%%%%%%%%%%%%%
+    % Load inputs if needed
+    %%%%%%%%%%%%%%%%%%%%%%%%%
+    
+    % Grain structure
+    if (~exist('grain', 'var') || (isempty(grain)))
+
+        disp('No grain structure found.');
+        phase_id = input('Specify the phase ID: \n', 's');
+        phase_id = str2double(phase_id);
+
+        grain_id = input('Specify the grain ID: \n', 's');
+        grain_id = str2double(grain_id);
+        grain = gtLoadGrain(phase_id, grain_id);
+    end
+
+    if (~exist('grain_rec', 'var') || (isempty(grain_rec)))
+
+        disp('No grain reconstruction details found.');
+        phase_id = input('Specify the phase ID: \n', 's');
+        phase_id = str2double(phase_id);
+
+        grain_id = input('Specify the grain ID: \n', 's');
+        grain_id = str2double(grain_id);
+        grain_rec = gtLoadGrainRec(phase_id, grain_id);
+
+    end
+
+    % Parameters
+    if (~exist('parameters', 'var') || isempty(parameters))
+
+        fprintf('No parameters file specified. Loading it from %s', pwd)
+        parameters = gtLoadParameters;
+
+    end
+    
+    % Deformation volume
+    if (~exist('dmvol', 'var') || isempty(dmvol))
+       
+       disp('No deformation volume specified.') 
+       phase_id = input('Specify the phase ID: \n', 's');
+       phase_id = str2double(phase_id);
+
+       grain_id = input('Specify the grain ID: \n', 's');
+       grain_id = str2double(grain_id);
+       grain_rec = gtLoadGrainRec(phase_id, grain_id);
+       dmvol = grain_rec.ODF6D.voxels_avg_R_vectors;
+
+    end
+
+    % Intensities volume
+    if (~exist('intvol', 'var') || isempty(intvol))
+
+       disp('No intensities volume specified.')
+       phase_id = input('Specify the phase ID: \n', 's');
+       phase_id = str2double(phase_id);
+
+       grain_id = input('Specify the grain ID: \n', 's');
+       grain_id = str2double(grain_id);
+       grain_rec = gtLoadGrainRec(phase_id, grain_id);
+       intvol = grain_rec.ODF6D.intensity;
+
+    end
+    
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    % Crop intensities and deformation volumes to save memory
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    disp('- Cropping intensities and deformation volumes ...')
+    tic
+    t = GtThreshold();
+    grain_rec
+    if (isfield(grain_rec, 'SEG'))
+        gr_seg = t.volumeThreshold(intvol, grain_rec.SEG.threshold, grain_rec.SEG.seed);
+    else
+        gr_seg = t.volumeThreshold(intvol, grain.threshold, grain.seed);
+    end
+    gr_seg
+    pause
+    cropbb = gr_seg.segbb + [1 1 1 0 0 0];
+    intvol = gtCrop(intvol, cropbb);
+
+    intvol(~gr_seg.seg) = 0;
+
+    dmvol_crop = zeros([size(intvol) 3]);
+    dmvol_crop(:, :, :, 1) = gtCrop(dmvol(:, :, :, 1), cropbb); 
+    dmvol_crop(:, :, :, 2) = gtCrop(dmvol(:, :, :, 2), cropbb);
+    dmvol_crop(:, :, :, 3) = gtCrop(dmvol(:, :, :, 3), cropbb);
+    dmvol = dmvol_crop;
+
+    toc()
+
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    % Creation of actual fedpars structure
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    disp('- Creating fedpars structure ...')
+
+    % x, y, z no. of voxels in the grain
+    fedpars.grainsize                  = size(intvol);
+
+    % x, y, z no. of voxels in a volume element
+    fedpars.gvsize                     = fedpars.grainsize ./ size(intvol); 
+%     gv_elem                        = numel(grain.vol);
+%     sample_freq                    = 3e+6;                       % verify
+%     gvsize                         = gv_elem / sample_freq;
+%     fedpars.gvsize                 = [gvsize gvsize gvsize]; 
+
+    % Active deformation components. Total number is 9.
+    % In case of rotations only (should be our case):
+    fedpars.dcomps                     = [true(1, 3) false(1, 6)];
+
+    % Origin voxel in the SAM reference
+    fedpars.vox000sam                  = parameters.samgeo.orig;
+
+    % Class of strain states
+    fedpars.gvtype                     = class(intvol);
+
+    % Real or simulated data?
+    fedpars.is_real_data               = true;
+    is_real = inputwdefault('Are you working with real data? [y/n]', 'y');
+
+    if (strcmpi(is_real, 'n'))
+        fedpars.is_real_data           = false;
+    end
+
+    % Generate random deformation components
+    if (~exist('dmvol', 'var') || isempty(dmvol))
+
+        % Type of deformation applied
+        fedpars.deftype                = 'smooth';
+        deftype = inputwdefault(['Specify the type of deformation ' ...
+            'applied: smooth (s) or realistic (r):', 's']);
+
+        if (strcomp(deftype, 'r'))
+            fedpars.deftype            = 'realistic';
+        end
+        
+        % Random strain definition mode
+        fedpars.randdefmode            = 'smooth';
+        randdefmode = input(['Specify the random strain definition mode: ' ...
+            'smooth (s) or interpolated (i):'], 's');
+
+        if (strcompi(randdefmode, 'i'))
+            fedpars.randdefmode        = 'interp';
+        end
+
+        % Maximum absolute deformation gradient
+        fedpars.dappgradlim            = 2e-4 * ones(size(fedpars.dcomps));    % guess: verify
+
+        % Maximum absolute deformation value which can be randomly generated
+        fedpars.dapplim                = 5e-3 * ones(size(fedpars.dcomps));    % guess: verify
+    end
+
+    % ID of blobs to load for figure
+    fedpars.loadbl = [];
+
+    % Detector parameters
+    ndet = numel(parameters.detgeo);
+
+    for ii_det = 1:ndet
+        
+        % Ondet
+        if (isfield(grain.proj, 'ondet') && (~isempty(grain.proj.ondet)))
+            fedpars.detector(ii_det).ondet = grain.proj.ondet;
+        else
+            fedpars.detector(ii_det).ondet = grain.ondet;
+        end
+
+        % Included
+        if (isfield(grain.proj, 'included') && (~isempty(grain.proj.included)))
+            fedpars.detector(ii_det).included = grain.proj.included;
+        else
+            fedpars.detector(ii_det).included = grain.included;
+        end
+
+        % Selected
+        if (isfield(grain.proj, 'selected') && (~isempty(grain.proj.selected)))
+            fedpars.detector(ii_det).selected = grain.proj.selected;
+        else
+            fedpars.detector(ii_det).selected = grain.selected;
+        end
+
+        % Increments in blobs bounding boxes accounting for deformation.
+        % Let gtFedTestLoadData_v2 decide the minimum size
+        fedpars.detector(ii_det).blobsizeadd = zeros(ndet, 3);
+
+    end
+
+    % Types for blobs intensities
+    fedpars.bltype                     = 'single';
+
+    % Deformation method for calculation of deformation tensor
+    fedpars.defmethod                  = 'small';
+
+    % Algorithm for blobs smoothing
+    fedpars.intsmoothalg               = 'none';
+
+    disp('Done')
+end
\ No newline at end of file

7_fed2/gtFedCreatePencilMask.m

+function intvol_out = gtFedCreatePencilMask(intvol, w, p, parameters)
+% intvol_out = gtPencilMask(intvol, w, p, parameters)
+%
+% Apply a pencil beam mask with square holes width x width (in microns) 
+% and pitch p (in microns) to an intensities volume reconstructed through
+% a 6D algorithm.
+
+    if (~exist('parameters', 'var') || isempty(parameters))
+        fprintf('No parameters file specified. Loading it from %s\n', pwd)
+        parameters = gtLoadParameters;
+    end
+
+    % Volume sizes in microns
+    xyz = size(intvol);
+    fprintf('Dimensions of intensities volume: %d %d %d\n', xyz(1), ...
+        xyz(2), xyz(3))
+
+    % Allocating space
+    mask = zeros(size(intvol));
+
+    w_vox = ceil((w / 1000) / parameters.recgeo.voxsize(2));
+    p_vox = ceil((p / 1000) / parameters.recgeo.voxsize(2));
+    fprintf(['%d central slices in the vertical direction are cut by '...
+        'the mask\n'], w_vox)
+
+    % Indices of vertical slices cut by the mask
+    v_ind(2) = floor((median(1:xyz(2))) + ceil(w_vox / 2));
+    v_ind(1) = v_ind(2) - w_vox + 1;
+
+    % Transmission function of the mask. The mask is periodic along 
+    % z-direction
+    cycle = [ones(1, w_vox) zeros(1, p_vox-w_vox)];
+    transmitted = repmat(cycle, [1, floor(xyz(3) / p_vox)]);
+    left = xyz(3) - numel(transmitted);
+    if (left <= w)
+        transmitted(numel(transmitted) + 1:xyz(3)) = ...
+               ones(1, xyz(3) - numel(transmitted));
+    else
+        transmitted(numel(transmitted) + 1:numel(transmitted) + w_vox) =...
+               ones(1, w_vox);
+        transmitted(numel(transmitted) + 1:xyz(3)) = ...
+               zeros(1, xyz(3) - numel(transmitted));   
+    end
+
+    transmitted = transmitted(ones(1, w_vox), :, ones(1, xyz(1)));
+    transmitted = permute(transmitted, [3 1 2]);
+    mask(:, v_ind(1):v_ind(2), :) = transmitted;
+    intvol_out = mask .* intvol;
+    fprintf('Dimensions of masked intensities volume: %d %d %d\n', ...
+        size(intvol_out, 1), size(intvol_out, 2), size(intvol_out, 3))
+end
\ No newline at end of file