From 3bc43db7cbdd87634836faeaa5df73ebb4edbd2d Mon Sep 17 00:00:00 2001
From: Lorenzo Valzania <valzania@rnice7-0308.esrf.fr>
Date: Wed, 4 Mar 2015 14:27:43 +0100
Subject: [PATCH] Creation of pencil mask for pencil beam illumination
Signed-off-by: Lorenzo Valzania <valzania@rnice7-0308.esrf.fr>
---
7_fed2/gtFedCreateGrainPars.m | 201 +++++++++++++++++++++++++++++++++
7_fed2/gtFedCreatePencilMask.m | 51 +++++++++
2 files changed, 252 insertions(+)
create mode 100644 7_fed2/gtFedCreateGrainPars.m
create mode 100644 7_fed2/gtFedCreatePencilMask.m
diff --git a/7_fed2/gtFedCreateGrainPars.m b/7_fed2/gtFedCreateGrainPars.m
new file mode 100644
index 00000000..868716c7
--- /dev/null
+++ b/7_fed2/gtFedCreateGrainPars.m
@@ -0,0 +1,201 @@
+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
diff --git a/7_fed2/gtFedCreatePencilMask.m b/7_fed2/gtFedCreatePencilMask.m
new file mode 100644
index 00000000..112c9d5b
--- /dev/null
+++ b/7_fed2/gtFedCreatePencilMask.m
@@ -0,0 +1,51 @@
+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
--
GitLab