From 8f3aecc8b52f1b4eb8f9c8618f9ff5fdb6cf5aa9 Mon Sep 17 00:00:00 2001 From: Nicola Vigano <nicola.vigano@esrf.fr> Date: Mon, 26 Feb 2018 18:35:49 +0100 Subject: [PATCH] 6D-reconstruction: got rid of useless dedicated path for super sampling It can be handled more easily the through the normal code path Signed-off-by: Nicola Vigano <nicola.vigano@esrf.fr> --- zUtil_Deformation/Gt6DBlobReconstructor.m | 28 +-- .../Gt6DReconstructionAlgorithmFactory.m | 175 +++++------------- zUtil_Deformation/Gt6DVolumeProjector.m | 64 +------ zUtil_Deformation/Gt6DVolumeToBlobProjector.m | 57 +----- 4 files changed, 66 insertions(+), 258 deletions(-) diff --git a/zUtil_Deformation/Gt6DBlobReconstructor.m b/zUtil_Deformation/Gt6DBlobReconstructor.m index 7003f613..29531425 100644 --- a/zUtil_Deformation/Gt6DBlobReconstructor.m +++ b/zUtil_Deformation/Gt6DBlobReconstructor.m @@ -345,18 +345,10 @@ classdef Gt6DBlobReconstructor < Gt6DVolumeToBlobProjector num_geoms = self.get_number_geometries(); num_det = self.get_number_detectors(); for ii_d = 1:num_det - if (self.using_super_sampling(ii_d)) - for n = 1:num_geoms - sinogram_size = [ ... - self.proj_size_uv(1), size(self.ss_geometries{ii_d}{n}, 1), self.proj_size_uv(2)]; - num_bytes = num_bytes + float_bytes * prod(sinogram_size); - end - else - for n = 1:num_geoms - sinogram_size = [ ... - self.proj_size_uv(1), size(self.geometries{ii_d}{n}, 1), self.proj_size_uv(2)]; - num_bytes = num_bytes + float_bytes * prod(sinogram_size); - end + for n = 1:num_geoms + sinogram_size = [ ... + self.proj_size_uv(1), size(self.geometries{ii_d}{n}, 1), self.proj_size_uv(2)]; + num_bytes = num_bytes + float_bytes * prod(sinogram_size); end end end @@ -692,16 +684,8 @@ classdef Gt6DBlobReconstructor < Gt6DVolumeToBlobProjector self.bwd_weights(1:num_ors) = {0}; for ii_d = 1:num_det - if (~self.using_super_sampling(ii_d)) - for ii = 1:num_ors - self.bwd_weights{ii} = self.bwd_weights{ii} + sum(self.offsets{ii_d}{ii}.proj_coeffs); - end - else - for ii = 1:num_ors - proj_structs = [self.ss_offsets{ii_d}{ii}{:}]; - proj_weights = cat(2, proj_structs(:).proj_coeffs); - self.bwd_weight{ii} = self.bwd_weights{ii} + sum(proj_weights); - end + for ii = 1:num_ors + self.bwd_weights{ii} = self.bwd_weights{ii} + sum(self.offsets{ii_d}{ii}.proj_coeffs); end end fprintf('\b\b (%2.1f s)\n', toc(c)); diff --git a/zUtil_Deformation/Gt6DReconstructionAlgorithmFactory.m b/zUtil_Deformation/Gt6DReconstructionAlgorithmFactory.m index 6db330e9..654398bb 100644 --- a/zUtil_Deformation/Gt6DReconstructionAlgorithmFactory.m +++ b/zUtil_Deformation/Gt6DReconstructionAlgorithmFactory.m @@ -70,8 +70,6 @@ classdef Gt6DReconstructionAlgorithmFactory < handle 'data_type', self.data_type, ... 'geometries', {algo_params.geometries}, ... 'offsets', {algo_params.offsets}, ... - 'ss_geometries', {algo_params.geometries_ss}, ... - 'ss_offsets', {algo_params.offsets_ss}, ... 'volume_ss', self.volume_downscaling, ... 'rspace_ss', self.rspace_super_sampling, ... 'orientation_groups', or_groups, ... @@ -186,10 +184,7 @@ classdef Gt6DReconstructionAlgorithmFactory < handle end geometries = cat(1, algo_params(:).geometries); - geometries_ss = cat(1, algo_params(:).geometries_ss); - offsets = algo_params(1).offsets; - offsets_ss = algo_params(1).offsets_ss; num_parent_blobs = numel(algo_params(1).blobs); @@ -262,22 +257,11 @@ classdef Gt6DReconstructionAlgorithmFactory < handle blobs{shared_blob_pos_inc} = new_blob; % Fixing the shift in parent's projection - if (~isempty(geometries_ss)) - for ii_o = 1:numel(offsets_ss) - for ii_ss = 1:numel(offsets_ss{ii_o}) - inds_of_shared_blobs = offsets_ss{ii_o}{ii_ss}.blob_offsets == shared_blob_pos_inc; - - offsets_ss{ii_o}{ii_ss}.proj_offsets(inds_of_shared_blobs) ... - = offsets_ss{ii_o}{ii_ss}.proj_offsets(inds_of_shared_blobs) + proj_shift; - end - end - else - for ii_o = 1:numel(offsets) - inds_of_shared_blobs = offsets{ii_o}.blob_offsets == shared_blob_pos_inc; + for ii_o = 1:numel(offsets) + inds_of_shared_blobs = offsets{ii_o}.blob_offsets == shared_blob_pos_inc; - offsets{ii_o}.proj_offsets(inds_of_shared_blobs) ... - = offsets{ii_o}.proj_offsets(inds_of_shared_blobs) + proj_shift; - end + offsets{ii_o}.proj_offsets(inds_of_shared_blobs) ... + = offsets{ii_o}.proj_offsets(inds_of_shared_blobs) + proj_shift; end end end @@ -298,79 +282,40 @@ classdef Gt6DReconstructionAlgorithmFactory < handle % Shifting projection to the concatenated blobs, unless that % blob was shared between parent and twin - if (~isempty(geometries_ss)) - proj_coeffs_ss = algo_params(ii_o_reg).offsets_ss; - for ii_o = 1:numel(proj_coeffs_ss) - for ii_ss = 1:numel(proj_coeffs_ss{ii_o}) - % Indices of the blobs, that we project to, in - % curr_or - inds_of_blobs = proj_coeffs_ss{ii_o}{ii_ss}.blob_offsets; - - % to_be_shifted contains the indices to the - % blobs belonging only to the current region - % identified by ii_o_reg - to_be_shifted = ~shared(inds_of_blobs); - - % The indices of the given blobs belonging to - % the current orientation will be shifted by - % the number of blobs accumulated already - inds_of_blobs(to_be_shifted) ... - = inds_of_blobs(to_be_shifted) + base_offset_shift; - - % Treating the shared blobs - inds_of_blobs(~to_be_shifted) ... - = ref_inc_to_sel(shared_pos(inds_of_blobs(~to_be_shifted))); - - proj_coeffs_ss{ii_o}{ii_ss}.blob_offsets = inds_of_blobs; - - % And now taking care of the positions inside - % those blobs - for ii_b = 1:numel(blobs_with_projs_to_be_shifted) - ind_of_parent_blob = parent_pos_shared_bls(blobs_with_projs_to_be_shifted(ii_b)); - inds_of_shared_blobs = inds_of_blobs == ind_of_parent_blob; - - proj_coeffs_ss{ii_o}{ii_ss}.proj_offsets(inds_of_shared_blobs) = ... - proj_coeffs_ss{ii_o}{ii_ss}.proj_offsets(inds_of_shared_blobs) + proj_shifts(blobs_with_projs_to_be_shifted(ii_b)); - end - end - end - offsets_ss = cat(1, offsets_ss, proj_coeffs_ss); - else - proj_coeffs = algo_params(ii_o_reg).offsets; - for ii_o = 1:numel(proj_coeffs) - % Indices of the blobs, that we project to, in - % curr_or - inds_of_blobs = proj_coeffs{ii_o}.blob_offsets; - - % to_be_shifted contains the indices to the - % blobs belonging only to the current region - % identified by ii_o_reg - to_be_shifted = ~shared(inds_of_blobs); - - % The indices of the given blobs belonging to - % the current orientation will be shifted by - % the number of blobs accumulated already - inds_of_blobs(to_be_shifted) ... - = inds_of_blobs(to_be_shifted) + base_offset_shift; - - % Treating the shared blobs - inds_of_blobs(~to_be_shifted) ... - = ref_inc_to_sel(shared_pos(inds_of_blobs(~to_be_shifted))); - - proj_coeffs{ii_o}.blob_offsets = inds_of_blobs; - - % And now taking care of the positions inside - % those blobs - for ii_b = 1:numel(blobs_with_projs_to_be_shifted) - ind_of_parent_blob = parent_pos_shared_bls(blobs_with_projs_to_be_shifted(ii_b)); - inds_of_shared_blobs = inds_of_blobs == ind_of_parent_blob; - - proj_coeffs{ii_o}.proj_offsets(inds_of_shared_blobs) = ... - proj_coeffs{ii_o}.proj_offsets(inds_of_shared_blobs) + proj_shifts(blobs_with_projs_to_be_shifted(ii_b)); - end + proj_coeffs = algo_params(ii_o_reg).offsets; + for ii_o = 1:numel(proj_coeffs) + % Indices of the blobs, that we project to, in + % curr_or + inds_of_blobs = proj_coeffs{ii_o}.blob_offsets; + + % to_be_shifted contains the indices to the + % blobs belonging only to the current region + % identified by ii_o_reg + to_be_shifted = ~shared(inds_of_blobs); + + % The indices of the given blobs belonging to + % the current orientation will be shifted by + % the number of blobs accumulated already + inds_of_blobs(to_be_shifted) ... + = inds_of_blobs(to_be_shifted) + base_offset_shift; + + % Treating the shared blobs + inds_of_blobs(~to_be_shifted) ... + = ref_inc_to_sel(shared_pos(inds_of_blobs(~to_be_shifted))); + + proj_coeffs{ii_o}.blob_offsets = inds_of_blobs; + + % And now taking care of the positions inside + % those blobs + for ii_b = 1:numel(blobs_with_projs_to_be_shifted) + ind_of_parent_blob = parent_pos_shared_bls(blobs_with_projs_to_be_shifted(ii_b)); + inds_of_shared_blobs = inds_of_blobs == ind_of_parent_blob; + + proj_coeffs{ii_o}.proj_offsets(inds_of_shared_blobs) = ... + proj_coeffs{ii_o}.proj_offsets(inds_of_shared_blobs) + proj_shifts(blobs_with_projs_to_be_shifted(ii_b)); end - offsets = cat(1, offsets, proj_coeffs); end + offsets = cat(1, offsets, proj_coeffs); base_offset_shift = base_offset_shift + numel(algo_params(ii_o_reg).blobs); end @@ -381,8 +326,6 @@ classdef Gt6DReconstructionAlgorithmFactory < handle 'data_type', self.data_type, ... 'geometries', geometries, ... 'offsets', offsets, ... - 'ss_geometries', geometries_ss, ... - 'ss_offsets', offsets_ss, ... 'volume_ss', self.volume_downscaling, ... 'psf', psf, ... 'orientation_groups', or_groups, ... @@ -606,17 +549,14 @@ classdef Gt6DReconstructionAlgorithmFactory < handle % Taking care of the orientation-space super-sampling tot_orient = sampler(ii_o_reg).get_total_num_orientations(); - if (num_ospace_oversampling == 1) - geometries_ss = {}; - offsets_ss = {}; - else + if (num_ospace_oversampling > 1) geometries = reshape(geometries, num_ospace_oversampling, []); offsets = reshape(offsets, num_ospace_oversampling, []); geometries_ss = cell(tot_orient, 1); offsets_ss = cell(tot_orient, 1); for ii_or = 1:tot_orient - geometries_ss{ii_or} = geometries(:, ii_or); + geometries_ss{ii_or} = cat(1, geometries{:, ii_or}); % merging all the orientation-space super-sampling % projection coefficients into one structure per @@ -625,21 +565,16 @@ classdef Gt6DReconstructionAlgorithmFactory < handle offsets_ss{ii_or} = self.merge_offsets(offsets(:, ii_or)); end - geometries = {}; - offsets = {}; + geometries = geometries_ss; + offsets = offsets_ss; end proj_props(ii_o_reg).geom = geometries; proj_props(ii_o_reg).offs = offsets; - proj_props(ii_o_reg).geom_ss = geometries_ss; - proj_props(ii_o_reg).offs_ss = offsets_ss; end geometries = cat(1, proj_props(:).geom); - geometries_ss = cat(1, proj_props(:).geom_ss); - offsets = cat(1, proj_props(:).offs); - offsets_ss = cat(1, proj_props(:).offs_ss); % Handling PSFs (very simplistic ATM) if (isfield(self.parameters.rec.grains.options, 'psf') ... @@ -659,19 +594,12 @@ classdef Gt6DReconstructionAlgorithmFactory < handle blobs_struct(1).det_ss = proj_det_ss; blobs_struct(1).psf = psf; - geometries = {geometries}; - geometries_ss = {geometries_ss}; - offsets = {offsets}; - offsets_ss = {offsets_ss}; - or_groups = self.get_orientation_groups(sampler); algo = Gt6DBlobReconstructor(volume_size, blobs_struct, ... 'data_type', self.data_type, ... - 'geometries', geometries, ... - 'offsets', offsets, ... - 'ss_geometries', geometries_ss, ... - 'ss_offsets', offsets_ss, ... + 'geometries', {geometries}, ... + 'offsets', {offsets}, ... 'volume_ss', self.volume_downscaling, ... 'orientation_groups', or_groups, ... varargin{:} ); @@ -859,17 +787,14 @@ classdef Gt6DReconstructionAlgorithmFactory < handle % Taking care of the orientation-space super-sampling tot_orient = sampler.get_total_num_orientations(); - if (num_ospace_oversampling == 1) - geometries_ss = {}; - offsets_ss = {}; - else + if (num_ospace_oversampling > 1) geometries = reshape(geometries, num_ospace_oversampling, []); offsets = reshape(offsets, num_ospace_oversampling, []); geometries_ss = cell(tot_orient, 1); offsets_ss = cell(tot_orient, 1); for ii_or = 1:tot_orient - geometries_ss{ii_or} = geometries(:, ii_or); + geometries_ss{ii_or} = cat(1, geometries{:, ii_or}); % merging all the orientation-space super-sampling % projection coefficients into one structure per @@ -878,8 +803,8 @@ classdef Gt6DReconstructionAlgorithmFactory < handle offsets_ss{ii_or} = self.merge_offsets(offsets(:, ii_or)); end - geometries = {}; - offsets = {}; + geometries = geometries_ss; + offsets = offsets_ss; end p = self.parameters; @@ -907,8 +832,6 @@ classdef Gt6DReconstructionAlgorithmFactory < handle 'blobs', {blobs_struct}, ... 'geometries', {geometries}, ... 'offsets', {offsets}, ... - 'geometries_ss', {geometries_ss}, ... - 'offsets_ss', {offsets_ss}, ... 'blobs_w_lims', {blobs_w_lims_pad} ); end @@ -1469,7 +1392,9 @@ classdef Gt6DReconstructionAlgorithmFactory < handle new_offsets_ss = offsets_ss{1}; cumulative_offs = max(offsets_ss{1}.sino_offsets); - for ii_ss = 2:numel(offsets_ss) + num_sub_orients = numel(offsets_ss); + + for ii_ss = 2:num_sub_orients new_offsets_ss.sino_offsets ... = cat(2, new_offsets_ss.sino_offsets, offsets_ss{ii_ss}.sino_offsets + cumulative_offs); new_offsets_ss.blob_offsets ... @@ -1482,6 +1407,8 @@ classdef Gt6DReconstructionAlgorithmFactory < handle num_projs = max(offsets_ss{ii_ss}.sino_offsets); cumulative_offs = cumulative_offs + num_projs; end + + new_offsets_ss.proj_coeffs = new_offsets_ss.proj_coeffs / num_sub_orients; end end end diff --git a/zUtil_Deformation/Gt6DVolumeProjector.m b/zUtil_Deformation/Gt6DVolumeProjector.m index 37a54deb..34a57461 100644 --- a/zUtil_Deformation/Gt6DVolumeProjector.m +++ b/zUtil_Deformation/Gt6DVolumeProjector.m @@ -10,8 +10,6 @@ classdef Gt6DVolumeProjector < handle astra_volume_geometry = {}; geometries = {}; - ss_geometries = {}; - astra_projection_geometries = {}; astra_projector_ids = {}; @@ -47,22 +45,11 @@ classdef Gt6DVolumeProjector < handle end num_det = self.get_number_detectors(); - if (isempty(self.ss_geometries{1})) - num_geoms = numel(self.geometries{1}); - else - num_geoms = numel(self.ss_geometries{1}); - end + num_geoms = numel(self.geometries{1}); for ii_d = 2:num_det - if (isempty(self.ss_geometries{1})) - if (num_geoms ~= numel(self.geometries{ii_d})) - error('Gt6DVolumeProjector:wrong_argument', ... - 'All detectors should have the same number of geometries!') - end - else - if (num_geoms ~= numel(self.ss_geometries{ii_d})) - error('Gt6DVolumeProjector:wrong_argument', ... - 'All detectors should have the same number of geometries!') - end + if (num_geoms ~= numel(self.geometries{ii_d})) + error('Gt6DVolumeProjector:wrong_argument', ... + 'All detectors should have the same number of geometries!') end end @@ -111,12 +98,7 @@ classdef Gt6DVolumeProjector < handle 'DetectorSuperSampling', self.rspace_ss, ... 'GPUindex', -1 ); - if (isempty(self.ss_geometries{1})) - num_geoms = numel(self.geometries{1}); - else - num_geoms = numel(self.ss_geometries{1}); - end - + num_geoms = numel(self.geometries{1}); num_det = self.get_number_detectors(); for ii_d = 1:num_det @@ -128,11 +110,7 @@ classdef Gt6DVolumeProjector < handle self.astra_projector_ids{ii_d} = cell(num_geoms, 1); for n = 1:num_geoms - if (isempty(self.ss_geometries{ii_d})) - geom = self.geometries{ii_d}{n}; - else - geom = cat(1, self.ss_geometries{ii_d}{n}{:}); - end + geom = self.geometries{ii_d}{n}; self.astra_projection_geometries{ii_d}{n} = astra_create_proj_geom(... 'parallel3d_vec', self.proj_size_uv(ii_d, 2), self.proj_size_uv(ii_d, 1), geom); @@ -151,10 +129,6 @@ classdef Gt6DVolumeProjector < handle num_geometries = numel(self.astra_projector_ids{1}); end - function is_ss = using_super_sampling(self, det_ind) - is_ss = ~isempty(self.ss_geometries{det_ind}); - end - function chunk_size = get_jobs_chunk_size(self) chunk_size = self.num_gpus * self.jobs_bunch_size; end @@ -172,10 +146,6 @@ classdef Gt6DVolumeProjector < handle function sinogram = fwd_project_volumes_to_sinos(self, volume, det_ind, n) % Basic Fwd-Projection function - if (self.using_super_sampling(det_ind)) - volume = self.rescale_volume_ss(volume, det_ind, n); - end - sinogram = astra_mex_direct_c('FP3D', [self.astra_projector_ids{det_ind}{n}], volume); end @@ -183,10 +153,6 @@ classdef Gt6DVolumeProjector < handle % Basic Bwd-Projection function volume = astra_mex_direct_c('BP3D', [self.astra_projector_ids{det_ind}{n}], sinogram); - - if (self.using_super_sampling(det_ind)) - volume = self.rescale_volume_ss(volume, det_ind, n); - end end function reset_geometry(self) @@ -203,22 +169,4 @@ classdef Gt6DVolumeProjector < handle end end end - - methods (Access = private) - function volume = rescale_single_volume_ss(self, volume, det_ind, n) - num_ss_geoms = numel(self.ss_geometries{det_ind}{n}); - volume = volume ./ num_ss_geoms; - end - - function volume = rescale_volume_ss(self, volume, det_ind, n) - num_orients = numel(n); - if (num_orients > 1) - for ii_o = 1:num_orients - volume{ii_o} = self.rescale_single_volume_ss(volume{ii_o}, det_ind, n(ii_o)); - end - else - volume = self.rescale_single_volume_ss(volume, det_ind, n); - end - end - end end diff --git a/zUtil_Deformation/Gt6DVolumeToBlobProjector.m b/zUtil_Deformation/Gt6DVolumeToBlobProjector.m index c533137f..9c3266ef 100644 --- a/zUtil_Deformation/Gt6DVolumeToBlobProjector.m +++ b/zUtil_Deformation/Gt6DVolumeToBlobProjector.m @@ -4,7 +4,6 @@ classdef Gt6DVolumeToBlobProjector < Gt6DVolumeProjector blobs; offsets = {}; - ss_offsets = {}; sinogram_sizes = zeros(0, 3, 0); @@ -41,11 +40,7 @@ classdef Gt6DVolumeToBlobProjector < Gt6DVolumeProjector num_det = self.get_number_detectors(); self.sinogram_sizes = zeros(num_orients, 3, num_det); for ii_d = 1:num_det - if (self.using_super_sampling(ii_d)) - geoms = self.ss_geometries{ii_d}; - else - geoms = self.geometries{ii_d}; - end + geoms = self.geometries{ii_d}; self.sinogram_sizes(:, :, ii_d) = [ ... self.proj_size_uv(ii_d * ones_n_ors, 1), zeros(num_orients, 1), self.proj_size_uv(ii_d * ones_n_ors, 2)]; @@ -54,44 +49,6 @@ classdef Gt6DVolumeToBlobProjector < Gt6DVolumeProjector end end end - -% function fwd_weights = getRowsSum(self) -% num_det = numel(self.blobs); -% -% chunk_size = self.get_jobs_chunk_size(); -% num_geoms = self.get_number_geometries(); -% -% vols_ones(1:chunk_size) = {ones(self.volume_geometry, self.data_type)}; -% -% fwd_weights = cell(num_det, 1); -% for ii_d = 1:num_det -% fwd_weights{ii_d} = gtMathsGetSameSizeZeros(self.blobs(ii_d).data); -% -% for n = 1:chunk_size:num_geoms -% chunk_safe_size = min(chunk_size, (num_geoms - n + 1)); -% -% inds = n:(n+chunk_safe_size-1); -% -% fwd_weights{ii_d} = self.fwd_project_volume( ... -% fwd_weights{ii_d}, vols_ones(1:chunk_safe_size), ii_d, inds); -% end -% end -% end - -% function bwd_weights = getColumnsSum(self) -% num_blobs = numel(self.blobs.sizes_w); -% blobs_ones = cell(1, num_blobs); -% for n = 1:num_blobs -% blobs_ones{n} = ones(self.proj_size_uv(1), ... -% self.blobs.sizes_w(n), self.proj_size_uv(2)); -% end -% -% num_geoms = self.get_number_geometries(); -% bwd_weights = cell(1, num_geoms); -% for n = 1:num_geoms -% bwd_weights{n} = self.bwd_project_volumes(blobs_ones, n); -% end -% end end methods (Access = protected) @@ -127,11 +84,7 @@ classdef Gt6DVolumeToBlobProjector < Gt6DVolumeProjector function [blobs, bs_time] = process_sinos_fwd(self, blobs, sinos, det_ind, ns) % Here we add the individual sinograms to the blobs - if (self.using_super_sampling(det_ind)) - offs = self.ss_offsets{det_ind}(ns); - else - offs = self.offsets{det_ind}(ns); - end + offs = self.offsets{det_ind}(ns); c = tic(); @@ -158,11 +111,7 @@ classdef Gt6DVolumeToBlobProjector < Gt6DVolumeProjector function [sinos, bs_time] = process_sinos_bwd(self, blobs, det_ind, ns) % Here we extract the individual sinograms from the blobs - if (self.using_super_sampling(det_ind)) - offs = self.ss_offsets{det_ind}(ns); - else - offs = self.offsets{det_ind}(ns); - end + offs = self.offsets{det_ind}(ns); sinos_size = self.sinogram_sizes(ns, :, det_ind); -- GitLab