diff --git a/zUtil_Deformation/Gt6DReconstructionAlgorithmFactory.m b/zUtil_Deformation/Gt6DReconstructionAlgorithmFactory.m index e594b99ffe073c29ec93cfc09226cbfdc910c032..11fc62de45d3a96c8c27d3460404ec80fbca2eb8 100644 --- a/zUtil_Deformation/Gt6DReconstructionAlgorithmFactory.m +++ b/zUtil_Deformation/Gt6DReconstructionAlgorithmFactory.m @@ -516,6 +516,8 @@ 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 @@ -524,42 +526,51 @@ classdef Gt6DReconstructionAlgorithmFactory < handle sf = sf_w{ii_o}; % Assign offsets - blob_offsets = cell(1, num_blobs); - proj_coeffs = cell(1, num_blobs); + ind_of_blob = cell(1, num_blobs); + slice_coeff = cell(1, num_blobs); ws = cell(1, num_blobs); for ii_b = 1:num_blobs - proj_coeffs{ii_b} = reshape(sf(ii_b).intm, 1, []); - - valid_coeffs = proj_coeffs{ii_b} > eps('single'); - proj_coeffs{ii_b} = proj_coeffs{ii_b}(valid_coeffs); + 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 - ws{ii_b} = sf(ii_b).bbwim(1):sf(ii_b).bbwim(2); - ws{ii_b} = ws{ii_b}(valid_coeffs); + valid_coeffs = cs > eps('single'); + slice_coeff{ii_b} = cs(valid_coeffs); + ws{ii_b} = iws(valid_coeffs); - blob_offsets{ii_b} = ii_b * ones(1, sum(valid_coeffs)); + ind_of_blob{ii_b} = ii_b * ones(1, sum(valid_coeffs)); end ws = [ws{:}]; - sino_offsets = 1:numel(ws); - blob_offsets = [blob_offsets{:}]; - proj_offsets = ws - lims_blobs_w(blob_offsets, 1)' + 1; - proj_coeffs = [proj_coeffs{:}]; + 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; + slice_coeff = [slice_coeff{:}]; offsets{ii_o} = struct( ... - 'sino_offsets', sino_offsets, ... - 'blob_offsets', blob_offsets, ... - 'proj_offsets', proj_offsets, ... - 'proj_coeffs', proj_coeffs ); + 'sino_offsets', pos_in_sino, ... + 'blob_offsets', ind_of_blob, ... + 'proj_offsets', pos_in_blob, ... + 'proj_coeffs', slice_coeff ); % Computing geometries rot_w_l2s = gtMathsRotationTensor(ws' .* om_step, rot_comp_w); - dvec_lab = ors_allblobs(ii_o).dvec(blob_offsets, :); + dvec_lab = ors_allblobs(ii_o).dvec(ind_of_blob, :); dvec_sam = gtGeoLab2Sam(dvec_lab, rot_w_l2s, ... self.parameters.labgeo, self.parameters.samgeo, true); proj_geom = gtGeoProjForReconstruction(... dvec_sam, rot_w_l2s, ref_gr.center, ... - bb_bls_uv(blob_offsets, :), [], ... + bb_bls_uv(ind_of_blob, :), [], ... self.parameters.detgeo(self.det_index), ... self.parameters.labgeo, ... self.parameters.samgeo, ... @@ -602,50 +613,82 @@ classdef Gt6DReconstructionAlgorithmFactory < handle bl_bbsizes = cat(1, sf.bbsize); % Assign offsets - blob_offsets = cell(1, num_blobs); - proj_coeffs = cell(1, num_blobs); + ind_of_blob = cell(1, num_blobs); + slice_coeff = cell(1, num_blobs); ns = cell(1, num_blobs); ws = cell(1, num_blobs); for ii_b = 1:num_blobs - ns{ii_b} = linspace(sf(ii_b).bbnim(1), sf(ii_b).bbnim(2), bl_bbsizes(ii_b, 1)); - ws{ii_b} = sf(ii_b).bbwim(1):sf(ii_b).bbwim(2); + 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 - ns{ii_b} = reshape(ns{ii_b}, [], 1); - ns{ii_b} = ns{ii_b}(:, ones(sf(ii_b).bbsize(2), 1)); - ns{ii_b} = reshape(ns{ii_b}, [], 1); + % 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)); - ws{ii_b} = ws{ii_b}(ones(sf(ii_b).bbsize(1), 1), :); - ws{ii_b} = reshape(ws{ii_b}, [], 1); + ins = reshape(ins, [], 1); + ins = ins(:, ones(numel(iws), 1)); + ins = reshape(ins, [], 1); - proj_coeffs{ii_b} = reshape(sf(ii_b).intm, 1, []); + iws = iws(ones(sf(ii_b).bbsize(1), 1), :); + iws = reshape(iws, [], 1); - valid_coeffs = proj_coeffs{ii_b} > eps('single'); - ns{ii_b} = ns{ii_b}(valid_coeffs); - ws{ii_b} = ws{ii_b}(valid_coeffs); - proj_coeffs{ii_b} = proj_coeffs{ii_b}(valid_coeffs); + % ins and iws have been expanded to be of the same + % size, but there is some redundancy going on there. + % This could impact on the perforance of the projection + % geometry determination later in this same function, + % but this is the only way to keep consistency with the + % sinogram<->blobs transformation coefficients. + valid_coeffs = cs > eps('single'); + ns{ii_b} = ins(valid_coeffs); + ws{ii_b} = iws(valid_coeffs); + slice_coeff{ii_b} = cs(valid_coeffs); - blob_offsets{ii_b} = ii_b * ones(1, sum(valid_coeffs)); + ind_of_blob{ii_b} = ii_b * ones(1, sum(valid_coeffs)); end ws = cat(1, ws{:}); ns = cat(1, ns{:}); - sino_offsets = 1:numel(ws); - blob_offsets = [blob_offsets{:}]; - proj_offsets = (ws - lims_blobs_w(blob_offsets, 1))' + 1; - proj_coeffs = [proj_coeffs{:}]; - - offsets{ii_o}.('sino_offsets') = sino_offsets; - offsets{ii_o}.('blob_offsets') = blob_offsets; - offsets{ii_o}.('proj_offsets') = proj_offsets; - offsets{ii_o}.('proj_coeffs') = proj_coeffs; + 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; + slice_coeff = [slice_coeff{:}]; - offsets{ii_o}.proj = struct( ... - 'sino_offset', sino_offsets, 'blob_offset', blob_offsets, ... - 'proj_offsets', proj_offsets, 'proj_coeffs', proj_coeffs ); + offsets{ii_o} = struct( ... + 'sino_offsets', pos_in_sino, ... + 'blob_offsets', ind_of_blob, ... + 'proj_offsets', pos_in_blob, ... + 'proj_coeffs', slice_coeff ); % Computing geometries rot_w_l2s = gtMathsRotationTensor(ws .* om_step, rot_comp_w); - t = ors_allblobs(ii_o).theta(blob_offsets); + % Building the projection geometry from the etas and omegas + t = ors_allblobs(ii_o).theta(ind_of_blob); pl_lab = gtGeoPlLabFromThetaEta(t, ns, self.parameters.labgeo); dvec_lab = gtFedPredictDiffVecMultiple(pl_lab', self.parameters.labgeo.beamdir')'; dvec_sam = gtGeoLab2Sam(dvec_lab, rot_w_l2s, ... @@ -653,7 +696,7 @@ classdef Gt6DReconstructionAlgorithmFactory < handle proj_geom = gtGeoProjForReconstruction(... dvec_sam, rot_w_l2s, ref_gr.center, ... - bb_bls_uv(blob_offsets, :), [], ... + bb_bls_uv(ind_of_blob, :), [], ... self.parameters.detgeo(self.det_index), ... self.parameters.labgeo, ... self.parameters.samgeo, ...