Skip to content
Snippets Groups Projects
Commit 65d2e378 authored by Nicola Vigano's avatar Nicola Vigano
Browse files

Shape-functions: reduced computational cost when using num_interp 


This also fixes a bug

Signed-off-by: default avatarNicola Vigano <vigano@yoda.esrf.fr>
parent 3d0bda7f
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
zUtil_Deformation/Gt6DReconstructionAlgorithmFactory.m
+ 8
− 43
View file @ 65d2e378
......@@ -705,7 +705,7 @@ classdef Gt6DReconstructionAlgorithmFactory < handle
end
function [geometries, offsets] = compute_proj_geom_shapefuncs_w(...
self, bl, ref_gr, ors_allblobs, lims_blobs_w, sf_w, slices_interp)
self, bl, ref_gr, ors_allblobs, lims_blobs_w, sf_w, blobs_w_interp)
% We here produce ASTRA's projection geometry and blob assempling
% coefficients
bls_bbuim = cat(1, bl.bbuim);
......@@ -720,8 +720,6 @@ classdef Gt6DReconstructionAlgorithmFactory < handle
om_step = gtGetOmegaStepDeg(self.parameters);
interp_w_blob_sizes = (lims_blobs_w(:, 2) - lims_blobs_w(:, 1)) ./ slices_interp + 1;
geometries = cell(num_orients, 1);
offsets = cell(num_orients, 1);
for ii_o = 1:num_orients
......@@ -735,17 +733,7 @@ classdef Gt6DReconstructionAlgorithmFactory < handle
ws = cell(1, num_blobs);
for ii_b = 1:num_blobs
cs = reshape(sf(ii_b).intm, 1, []);
iws = sf(ii_b).bbwim(1):sf(ii_b).bbwim(2);
if (slices_interp(ii_b) > 1)
ni_pos = (iws - lims_blobs_w(ii_b, 1)) ./ slices_interp(ii_b) + 1;
[inds, ints] = gtMathsGetInterpolationIndices(ni_pos', cs');
valid_inds = ints > eps('single');
inds = inds(valid_inds);
ints = ints(valid_inds);
cs = accumarray(inds, ints, [interp_w_blob_sizes(ii_b), 1])';
iws = (0:numel(cs)-1) .* slices_interp(ii_b) + lims_blobs_w(ii_b, 1);
end
iws = sf(ii_b).bbwim(1):blobs_w_interp(ii_b):sf(ii_b).bbwim(2);
valid_coeffs = cs > eps('single');
slice_coeff{ii_b} = cs(valid_coeffs);
......@@ -753,10 +741,11 @@ classdef Gt6DReconstructionAlgorithmFactory < handle
ind_of_blob{ii_b} = ii_b * ones(1, sum(valid_coeffs));
end
ws = [ws{:}];
pos_in_sino = 1:numel(ws);
ind_of_blob = [ind_of_blob{:}];
pos_in_blob = (ws - lims_blobs_w(ind_of_blob, 1)') ./ slices_interp(ind_of_blob, 1)' + 1;
pos_in_blob = (ws - lims_blobs_w(ind_of_blob, 1)') ./ blobs_w_interp(ind_of_blob, 1)' + 1;
slice_coeff = [slice_coeff{:}];
offsets{ii_o} = struct( ...
......@@ -792,7 +781,7 @@ classdef Gt6DReconstructionAlgorithmFactory < handle
end
function [geometries, offsets] = compute_proj_geom_shapefuncs_nw(...
self, bl, ref_gr, ors_allblobs, lims_blobs_w, sf_nw, slices_interp)
self, bl, ref_gr, ors_allblobs, lims_blobs_w, sf_nw, blobs_w_interp)
% We here produce ASTRA's projection geometry and blob assempling
% coefficients
bls_bbuim = cat(1, bl.bbuim);
......@@ -825,32 +814,7 @@ classdef Gt6DReconstructionAlgorithmFactory < handle
cs = sf(ii_b).intm;
cs = reshape(cs, 1, []);
% Omegas of each column of the shape-function
iws = sf(ii_b).bbwim(1):sf(ii_b).bbwim(2);
if (slices_interp(ii_b) > 1)
% Compressing the Omegas, to fit in num_interp
iws_numint = (iws - lims_blobs_w(ii_b, 1)) ./ slices_interp(ii_b) + 1;
iws_numint = iws_numint(ones(sf(ii_b).bbsize(1), 1), :);
iws_numint = reshape(iws_numint, [], 1);
ins_numint = reshape(1:sf(ii_b).bbsize(1), [], 1);
ins_numint = ins_numint(:, ones(numel(iws), 1));
ins_numint = reshape(ins_numint, [], 1);
[inds, ints] = gtMathsGetInterpolationIndices([ins_numint, iws_numint], cs');
inds = inds(:, 1) + (inds(:, 2) - 1) * sf(ii_b).bbsize(1);
valid_inds = ints > eps('single');
inds = inds(valid_inds);
ints = ints(valid_inds);
% We now create the new coefficients, having some 0
% entries in the lower side of the w directions,
% which will then filtered out by the valid_coeffs
% variable.
size_w_sf = max(ceil(iws_numint));
cs = accumarray(inds, ints, [sf(ii_b).bbsize(1) * size_w_sf, 1])';
iws = (0:size_w_sf-1) .* slices_interp(ii_b) + lims_blobs_w(ii_b, 1);
end
iws = sf(ii_b).bbwim(1):blobs_w_interp(ii_b):sf(ii_b).bbwim(2);
% Etas of each row of the shape-function
ins = linspace(sf(ii_b).bbnim(1), sf(ii_b).bbnim(2), bl_bbsizes(ii_b, 1));
......@@ -875,11 +839,12 @@ classdef Gt6DReconstructionAlgorithmFactory < handle
ind_of_blob{ii_b} = ii_b * ones(1, sum(valid_coeffs));
end
ws = cat(1, ws{:});
ns = cat(1, ns{:});
pos_in_sino = 1:numel(ws);
ind_of_blob = [ind_of_blob{:}];
pos_in_blob = ((ws - lims_blobs_w(ind_of_blob, 1)) ./ slices_interp(ind_of_blob, 1))' + 1;
pos_in_blob = ((ws - lims_blobs_w(ind_of_blob, 1)) ./ blobs_w_interp(ind_of_blob, 1))' + 1;
slice_coeff = [slice_coeff{:}];
offsets{ii_o} = struct( ...
......
zUtil_Deformation/gtDefFwdProjGvdm2NW.m
+ 21
− 6
View file @ 65d2e378
......@@ -33,10 +33,23 @@ function bl = gtDefFwdProjGvdm2NW(grain, ref_sel, gv, fedpars, parameters, eta_s
pls_orig = grain.allblobs.plorig(ref_sel, :);
end
if (isfield(fedpars, 'detector') ...
&& isfield(fedpars.detector(det_ind), 'blobs_w_interp') ...
&& ~isempty(fedpars.detector(det_ind).blobs_w_interp))
blobs_w_interp = fedpars.detector(det_ind).blobs_w_interp;
if (numel(blobs_w_interp) ~= nbl)
error('gtDefFwdProjGvdm2W:wrong_argument', ...
'Number of "blobs_w_interp" (%d) doesn''t match with the number of blobs (%d)', ...
numel(blobs_w_interp), nbl)
end
else
blobs_w_interp = ones(nbl, 1);
end
sinths = grain.allblobs.sintheta(ref_sel);
ominds = grain.allblobs.omind(ref_sel);
ref_ws = grain.allblobs.detector(det_ind).uvw(ref_sel, 3);
w_shifts = round(ref_ws);
w_shifts = round(ref_ws ./ blobs_w_interp);
ref_ns = grain.allblobs.eta(ref_sel);
n_shifts = round(ref_ns ./ eta_step);
......@@ -95,15 +108,16 @@ function bl = gtDefFwdProjGvdm2NW(grain, ref_sel, gv, fedpars, parameters, eta_s
numel(inds_bad), numel(om), ii_b)
end
nw_bl = om ./ om_step;
nw_bl = gtGrainAnglesTabularFix360deg(nw_bl, ref_ws(ii_b), parameters);
nw_bl = nw_bl - w_shifts(ii_b);
w_bl = om ./ om_step;
w_bl = gtGrainAnglesTabularFix360deg(w_bl, ref_ws(ii_b), parameters);
w_bl = w_bl / blobs_w_interp(ii_b);
w_bl = w_bl - w_shifts(ii_b);
n_bl = gtGeoEtaFromDiffVec(pl_lab', parameters.labgeo)';
n_bl = gtGrainAnglesTabularFix360deg(n_bl, ref_ns(ii_b));
n_bl = n_bl ./ eta_step - n_shifts(ii_b);
nw{ii_b} = [n_bl; nw_bl];
nw{ii_b} = [n_bl; w_bl];
nw_min(ii_b, :) = min(nw{ii_b}, [], 2);
nw_max(ii_b, :) = max(nw{ii_b}, [], 2);
......@@ -181,7 +195,8 @@ function bl = gtDefFwdProjGvdm2NW(grain, ref_sel, gv, fedpars, parameters, eta_s
bl(ii_b).bbnim = ([nw_min(ii_b, 1), nw_max(ii_b, 1)] + n_shifts(ii_b)) .* eta_step;
im_w_low_lim = w_shifts(ii_b) - blob_orig_nw_shifts(ii_b, 2) + 1;
bl(ii_b).bbwim = [im_w_low_lim, im_w_low_lim + bl(ii_b).bbsize(2) - 1];
bl(ii_b).bbwim = [im_w_low_lim, im_w_low_lim + bl(ii_b).bbsize(2) - 1] * blobs_w_interp(ii_b);
bl(ii_b).interpw = blobs_w_interp(ii_b);
o.fprintf(repmat('\b', [1, num_chars]));
end
......
zUtil_Deformation/gtDefFwdProjGvdm2W.m
+ 17
− 2
View file @ 65d2e378
......@@ -33,10 +33,23 @@ function bl = gtDefFwdProjGvdm2W(grain, ref_sel, gv, fedpars, parameters, det_in
pls_orig = grain.allblobs.plorig(ref_sel, :);
end
if (isfield(fedpars, 'detector') ...
&& isfield(fedpars.detector(det_ind), 'blobs_w_interp') ...
&& ~isempty(fedpars.detector(det_ind).blobs_w_interp))
blobs_w_interp = fedpars.detector(det_ind).blobs_w_interp;
if (numel(blobs_w_interp) ~= nbl)
error('gtDefFwdProjGvdm2W:wrong_argument', ...
'Number of "blobs_w_interp" (%d) doesn''t match with the number of blobs (%d)', ...
numel(blobs_w_interp), nbl)
end
else
blobs_w_interp = ones(nbl, 1);
end
sinths = grain.allblobs.sintheta(ref_sel);
ominds = grain.allblobs.omind(ref_sel);
ref_ws = grain.allblobs.detector(det_ind).uvw(ref_sel, 3);
w_shifts = round(ref_ws);
w_shifts = round(ref_ws ./ blobs_w_interp);
% The plane normals need to be brought in the Lab reference where the
% beam direction and rotation axis are defined.
......@@ -95,6 +108,7 @@ function bl = gtDefFwdProjGvdm2W(grain, ref_sel, gv, fedpars, parameters, det_in
w_bl = om ./ om_step;
w_bl = gtGrainAnglesTabularFix360deg(w_bl, ref_ws(ii_b), parameters);
w_bl = w_bl / blobs_w_interp(ii_b);
w_bl = w_bl - w_shifts(ii_b);
w{ii_b} = w_bl;
......@@ -173,7 +187,8 @@ function bl = gtDefFwdProjGvdm2W(grain, ref_sel, gv, fedpars, parameters, det_in
bl(ii_b).bbsize = blob_w_sizes(ii_b);
im_low_lim = w_shifts(ii_b) - blob_orig_w_shifts(ii_b) + 1;
bl(ii_b).bbwim = [im_low_lim, im_low_lim + bl(ii_b).bbsize - 1];
bl(ii_b).bbwim = [im_low_lim, im_low_lim + bl(ii_b).bbsize - 1] * blobs_w_interp(ii_b);
bl(ii_b).interpw = blobs_w_interp(ii_b);
o.fprintf(repmat('\b', [1, num_chars]));
end
......
zUtil_Deformation/gtDefShapeFunctionsFwdProj.m
+ 17
− 3
View file @ 65d2e378
......@@ -10,7 +10,7 @@ function shape_funcs = gtDefShapeFunctionsFwdProj(sampler, varargin)
'interpolated_voxel', false, ...
'num_interp', [], ...
'blobs_w_interp', [], ...
'projector', 'nearest');
'projector', 'linear');
conf = parse_pv_pairs(conf, varargin);
num_detectors = numel(sampler.parameters.detgeo);
......@@ -20,8 +20,18 @@ function shape_funcs = gtDefShapeFunctionsFwdProj(sampler, varargin)
labgeo.rotcomp = gtMathsRotationMatrixComp(labgeo.rotdir', 'col');
sampler.parameters.labgeo = labgeo;
use_numinterp = ~isempty(conf.num_interp) && ~isempty(conf.blobs_w_interp);
fed_pars_detector = struct(...
'blobsizeadd', [0 0 0], 'psf', [], 'apply_uv_shift', conf.recenter_sf);
'blobsizeadd', [0 0 0], ...
'psf', [], ...
'apply_uv_shift', conf.recenter_sf, ...
'blobs_w_interp', []);
if (use_numinterp)
fed_pars_detector.blobs_w_interp = conf.blobs_w_interp;
end
fedpars = struct( ...
'bltype', conf.data_type, ...
'dcomps', [1 1 1 0 0 0 0 0 0] == 1, ...
......@@ -32,7 +42,11 @@ function shape_funcs = gtDefShapeFunctionsFwdProj(sampler, varargin)
voxel_size = sampler.stats.sampling.gaps;
space_res = tand(sampler.estimate_maximum_resolution() ./ 2);
voxel_sampling = max(ceil(voxel_size ./ space_res * conf.factor), 5);
voxel_sampling = voxel_size ./ space_res * conf.factor;
if (use_numinterp)
voxel_sampling(3) = voxel_sampling(3) ./ conf.num_interp;
end
voxel_sampling = max(ceil(voxel_sampling), 5);
num_sub_voxels = prod(voxel_sampling);
ones_nsv = ones(num_sub_voxels, 1);
......
zUtil_ForwardSim/gtFwdSimBlobDefinition.m
+ 4
− 2
View file @ 65d2e378
......@@ -11,7 +11,8 @@ function blob = gtFwdSimBlobDefinition(type, num_blobs)
blob = struct( ...
'intm', [], ... % The normalized shape function
'bbwim', zeros(0, 2), ... % Shape function's limits in W
'bbsize', zeros(0, 1) ); % Shape function's BoundingBox (includes margins)
'bbsize', zeros(0, 1), ... % Shape function's BoundingBox (includes margins)
'interpw', zeros(0, 1) ); % Size of the W interpolation (num_interp)
case 'sf_n'
blob = struct( ...
'intm', [], ... % The normalized shape function
......@@ -22,7 +23,8 @@ function blob = gtFwdSimBlobDefinition(type, num_blobs)
'intm', [], ... % The normalized shape function
'bbnim', zeros(0, 2), ... % Shape function's limits in eta
'bbwim', zeros(0, 2), ... % Shape function's limits in W
'bbsize', zeros(0, 2) ); % Shape function's BoundingBox (includes margins)
'bbsize', zeros(0, 2), ... % Shape function's BoundingBox (includes margins)
'interpw', zeros(0, 1) ); % Size of the W interpolation (num_interp)
case {'blob', 'sf_uvw'}
blob = struct( ...
'intm', [], ... % The normalized blob
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment