diff --git a/zUtil_Deformation/Gt6DBlobReconstructor.m b/zUtil_Deformation/Gt6DBlobReconstructor.m index 87a2ff62d2269acd765f457a4da913c39812f100..12ab76c1818cdbc3face82cba126059f08b31ab4 100644 --- a/zUtil_Deformation/Gt6DBlobReconstructor.m +++ b/zUtil_Deformation/Gt6DBlobReconstructor.m @@ -39,7 +39,6 @@ classdef Gt6DBlobReconstructor < Gt6DVolumeToBlobProjector self.statistics.add_task('cp_dual_update_detector', 'CP Dual variable (detector) update'); self.statistics.add_task_partial('cp_dual_update_detector', 'cp_dual_detector_FP', 'Forward Projection'); self.statistics.add_task_partial('cp_dual_update_detector', 'cp_dual_detector_SB', 'Sinograms -> Blobs'); - self.statistics.add_task_partial('cp_dual_update_detector', 'cp_dual_detector_SF', 'Shape Functions'); self.statistics.add_task('cp_dual_update_l1', 'CP Dual variable (l1) update'); self.statistics.add_task('cp_dual_update_tv', 'CP Dual variable (TV) update'); self.statistics.add_task_partial('cp_dual_update_tv', 'cp_dual_tv_reduction', 'Volumes reduction'); @@ -50,7 +49,6 @@ classdef Gt6DBlobReconstructor < Gt6DVolumeToBlobProjector self.statistics.add_task('cp_primal_update', 'CP Primal variable update'); self.statistics.add_task_partial('cp_primal_update', 'cp_primal_BP', 'Back Projection'); self.statistics.add_task_partial('cp_primal_update', 'cp_primal_BS', 'Blobs -> Sinograms'); - self.statistics.add_task_partial('cp_primal_update', 'cp_primal_SF', 'Transposed Shape Functions'); self.statistics.add_task_partial('cp_primal_update', 'cp_primal_ODF', 'ODF Correction'); self.statistics.add_task_partial('cp_primal_update', 'cp_primal_CORR', 'Primal correction computation'); self.statistics.add_task_partial('cp_primal_update', 'cp_primal_APP', 'Primal update application'); @@ -410,7 +408,6 @@ classdef Gt6DBlobReconstructor < Gt6DVolumeToBlobProjector timing_fp = 0; timing_sb = 0; - timing_sf_fp = 0; num_geoms = self.get_number_geometries(); @@ -421,30 +418,27 @@ classdef Gt6DBlobReconstructor < Gt6DVolumeToBlobProjector for n = 1:chunk_size:chunk_end inds = n:(n+chunk_size-1); - [proj_bls, fp_time, sb_time, sf_fp_time] = self.fwd_project_volume( ... + [proj_bls, fp_time, sb_time] = self.fwd_project_volume( ... proj_bls, volumes(inds), inds); timing_fp = timing_fp + fp_time; timing_sb = timing_sb + sb_time; - timing_sf_fp = timing_sf_fp + sf_fp_time; end else chunk_end = 0; end for n = (chunk_end+1):num_geoms - [proj_bls, fp_time, sb_time, sf_fp_time] = self.fwd_project_volume( ... + [proj_bls, fp_time, sb_time] = self.fwd_project_volume( ... proj_bls, volumes{n}, n); timing_fp = timing_fp + fp_time; timing_sb = timing_sb + sb_time; - timing_sf_fp = timing_sf_fp + sf_fp_time; end if (apply_stats) self.statistics.add_timestamp(timing_fp, 'cp_dual_update_detector', 'cp_dual_detector_FP'); self.statistics.add_timestamp(timing_sb, 'cp_dual_update_detector', 'cp_dual_detector_SB'); - self.statistics.add_timestamp(timing_sf_fp, 'cp_dual_update_detector', 'cp_dual_detector_SF'); end end @@ -453,7 +447,6 @@ classdef Gt6DBlobReconstructor < Gt6DVolumeToBlobProjector timing_bs = 0; timing_corr = 0; timing_app = 0; - timing_sf_bp = 0; num_geoms = self.get_number_geometries(); @@ -464,10 +457,9 @@ classdef Gt6DBlobReconstructor < Gt6DVolumeToBlobProjector for n = 1:chunk_size:chunk_end inds = n:(n+chunk_size-1); - [v, bp_time, bs_time, sf_bp_time] = self.bwd_project_volume(p_blurr, inds); + [v, bp_time, bs_time] = self.bwd_project_volume(p_blurr, inds); timing_bp = timing_bp + bp_time; timing_bs = timing_bs + bs_time; - timing_sf_bp = timing_sf_bp + sf_bp_time; for ii_gpu = 1:chunk_size % Computing the update to apply @@ -486,10 +478,9 @@ classdef Gt6DBlobReconstructor < Gt6DVolumeToBlobProjector end for n = (chunk_end+1):num_geoms - [v, bp_time, bs_time, sf_bp_time] = self.bwd_project_volume(p_blurr, n); + [v, bp_time, bs_time] = self.bwd_project_volume(p_blurr, n); timing_bp = timing_bp + bp_time; timing_bs = timing_bs + bs_time; - timing_sf_bp = timing_sf_bp + sf_bp_time; % Computing the update to apply up_prim = [compute_update_primal(n), v]; @@ -503,7 +494,6 @@ classdef Gt6DBlobReconstructor < Gt6DVolumeToBlobProjector self.statistics.add_timestamp(timing_bp, 'cp_primal_update', 'cp_primal_BP') self.statistics.add_timestamp(timing_bs, 'cp_primal_update', 'cp_primal_BS') - self.statistics.add_timestamp(timing_sf_bp, 'cp_primal_update', 'cp_primal_SF') self.statistics.add_timestamp(timing_corr, 'cp_primal_update', 'cp_primal_CORR') self.statistics.add_timestamp(timing_app, 'cp_primal_update', 'cp_primal_APP') end diff --git a/zUtil_Deformation/Gt6DReconstructionAlgorithmFactory.m b/zUtil_Deformation/Gt6DReconstructionAlgorithmFactory.m index 17494f8ad4c530dc706fcb37417d2b0d37c94b66..c09dc540be8b9bdf8cd533f0ce24fab22768375c 100644 --- a/zUtil_Deformation/Gt6DReconstructionAlgorithmFactory.m +++ b/zUtil_Deformation/Gt6DReconstructionAlgorithmFactory.m @@ -16,8 +16,6 @@ classdef Gt6DReconstructionAlgorithmFactory < handle % 'none' : no shape functions are used % 'w' : Only omega is considered % 'nw' : (eta, omega) shape functions - % 'nw2uvw' : (eta, omega) shape functions -> UVW shape functions - % 'uvw' : Using FED's projector for UVW shape functions shape_functions_type = 'none'; parameters; @@ -47,23 +45,13 @@ classdef Gt6DReconstructionAlgorithmFactory < handle switch (lower(self.shape_functions_type)) case 'none' - shape_funcs = {}; algo_params = self.get_algorithm_params_num_interp(sampler, num_interp); case 'w' shape_funcs = gtDefShapeFunctionsCreateW(sampler); algo_params = self.get_algorithm_params_shape_funcs(sampler, shape_funcs, self.shape_functions_type); - shape_funcs = {}; case 'nw' shape_funcs = gtDefShapeFunctionsCreateNW(sampler); algo_params = self.get_algorithm_params_shape_funcs(sampler, shape_funcs, self.shape_functions_type); - shape_funcs = {}; - case 'nw2uvw' - shape_funcs = gtDefShapeFunctionsCreateNW(sampler, 2, 'data_type', 'single'); - shape_funcs = gtDefShapeFunctionsNW2UVW(sampler, shape_funcs); - algo_params = self.get_algorithm_params_shape_funcs(sampler, shape_funcs, self.shape_functions_type); - case 'uvw' - shape_funcs = gtDefShapeFunctionsFwdProjUVW(sampler, 'factor', 15); - algo_params = self.get_algorithm_params_shape_funcs(sampler, shape_funcs, self.shape_functions_type); otherwise error('Gt6DReconstructionAlgorithmFactory:wrong_argument', ... 'Shape functions of type: %s, not supported yet!', ... @@ -82,7 +70,6 @@ classdef Gt6DReconstructionAlgorithmFactory < handle 'volume_ss', self.volume_downscaling, ... 'rspace_ss', self.rspace_super_sampling, ... 'psf', algo_params.psf, ... - 'shape_functions', shape_funcs, ... varargin{:} ); end @@ -466,10 +453,6 @@ classdef Gt6DReconstructionAlgorithmFactory < handle case 'nw' [geometries, offsets] = self.compute_proj_geom_shapefuncs_nw( ... bl, ref_gr, ors_allblobs, shape_funcs, lims_blobs_w); - case {'nw2uvw', 'uvw'} - geometries = self.compute_proj_geom_numinterp(bl, ref_gr, ors_allblobs); - offsets = self.make_proj_geom_shape_funcs_uvw(... - shape_funcs, lims_blobs_w, w_tab); end tot_orient = sampler.get_total_num_orientations(); @@ -1030,47 +1013,6 @@ classdef Gt6DReconstructionAlgorithmFactory < handle end end - function offsets = make_proj_geom_shape_funcs_uvw(self, shape_funcs, lims_blobs_w, w_tab) - num_geoms = size(w_tab, 2); - num_blobs = size(w_tab, 1); - - slice_pos_in_blob = w_tab - repmat(lims_blobs_w(:, 1), [1 num_geoms 2]) + 1; - - % Assign offsets - offsets = cell(num_geoms, 1); - for ii = 1:num_geoms - sino_offsets = cell(1, num_blobs); - blob_offsets = cell(1, num_blobs); - proj_offsets = cell(1, num_blobs); - proj_coeffs = cell(1, num_blobs); - - cumul_pos = 0; - - for n = 1:num_blobs - shape_func = shape_funcs{ii}(n); - ones_num_slices = ones(1, shape_func.bbsize(3)); - - slices = slice_pos_in_blob(n, ii, 1):slice_pos_in_blob(n, ii, 2); - sinos_offs = cumul_pos + (1:shape_func.bbsize(3)); - cumul_pos = cumul_pos + shape_func.bbsize(3); - - sino_offsets{n} = sinos_offs; - blob_offsets{n} = n * ones_num_slices; - proj_offsets{n} = slices; - proj_coeffs{n} = ones_num_slices; - end - - offsets{ii}.('sino_offsets') = [sino_offsets{:}]; - offsets{ii}.('blob_offsets') = [blob_offsets{:}]; - offsets{ii}.('proj_offsets') = [proj_offsets{:}]; - offsets{ii}.('proj_coeffs') = [proj_coeffs{:}]; - - offsets{ii}.proj = struct( ... - 'sino_offset', sino_offsets, 'blob_offset', blob_offsets, ... - 'proj_offsets', proj_offsets, 'proj_coeffs', proj_coeffs ); - end - end - function blobs = renormalize_blobs(self, sampler, blobs) ref_gr = sampler.get_reference_grain(); if (isfield(self.parameters.cryst(ref_gr.phaseid), 'int_exp')) diff --git a/zUtil_Deformation/Gt6DVolumeToBlobProjector.m b/zUtil_Deformation/Gt6DVolumeToBlobProjector.m index c942ce4091567cd5ba315a4d25b89535b0ecd9a8..da1a52ed516805a83a206f64f50984bbab8239a0 100644 --- a/zUtil_Deformation/Gt6DVolumeToBlobProjector.m +++ b/zUtil_Deformation/Gt6DVolumeToBlobProjector.m @@ -8,8 +8,6 @@ classdef Gt6DVolumeToBlobProjector < Gt6DVolumeProjector fwd_weights = {}; bwd_weights = {}; - shape_functions = {}; - blob_sinogram_c_functions = true; end @@ -68,7 +66,7 @@ classdef Gt6DVolumeToBlobProjector < Gt6DVolumeProjector end methods (Access = protected) - function [blobs, fp_time, bs_time, sf_time] = fwd_project_volume(self, blobs, volume, n) + function [blobs, fp_time, bs_time] = fwd_project_volume(self, blobs, volume, n) % We get volumes in input and we want the blobs as output, but the % output from the base function is stacks of sinograms relative to % the geometry of the volumes, so we need to actually decompose the @@ -81,18 +79,16 @@ classdef Gt6DVolumeToBlobProjector < Gt6DVolumeProjector num_orients = numel(n); if (num_orients > 1) bs_time = 0; - sf_time = 0; for ii_o = 1:num_orients - [blobs, bs_time_o, sf_time_o] = self.process_sino_fwd(blobs, sino{ii_o}, n(ii_o)); + [blobs, bs_time_o] = self.process_sino_fwd(blobs, sino{ii_o}, n(ii_o)); bs_time = bs_time + bs_time_o; - sf_time = sf_time + sf_time_o; end else - [blobs, bs_time, sf_time] = self.process_sino_fwd(blobs, sino, n); + [blobs, bs_time] = self.process_sino_fwd(blobs, sino, n); end end - function [volume, bp_time, bs_time, sf_time] = bwd_project_volume(self, blobs, n) + function [volume, bp_time, bs_time] = bwd_project_volume(self, blobs, n) % Similar problem to the projection function, with the only % difference that in this case the input is blobs and we want to % transform them into stack of sinograms @@ -101,14 +97,12 @@ classdef Gt6DVolumeToBlobProjector < Gt6DVolumeProjector if (num_orients > 1) sino = cell(num_orients, 1); bs_time = 0; - sf_time = 0; for ii_o = 1:num_orients - [sino{ii_o}, bs_time_o, sf_time_o] = self.process_sino_bwd(blobs, n(ii_o)); + [sino{ii_o}, bs_time_o] = self.process_sino_bwd(blobs, n(ii_o)); bs_time = bs_time + bs_time_o; - sf_time = sf_time + sf_time_o; end else - [sino, bs_time, sf_time] = self.process_sino_bwd(blobs, n); + [sino, bs_time] = self.process_sino_bwd(blobs, n); end c = tic(); @@ -116,17 +110,9 @@ classdef Gt6DVolumeToBlobProjector < Gt6DVolumeProjector bp_time = toc(c); end - function [blobs, bs_time, sf_time] = process_sino_fwd(self, blobs, sino, n) + function [blobs, bs_time] = process_sino_fwd(self, blobs, sino, n) % Here we add the individual sinograms to the blobs - if (~isempty(self.shape_functions)) - c = tic(); - sino = self.apply_shape_functions_fp(sino, self.shape_functions{n}); - sf_time = toc(c); - else - sf_time = []; - end - if (self.using_super_sampling()) offs = self.ss_offsets{n}; else @@ -144,7 +130,7 @@ classdef Gt6DVolumeToBlobProjector < Gt6DVolumeProjector bs_time = toc(c); end - function [sino, bs_time, sf_time] = process_sino_bwd(self, blobs, n) + function [sino, bs_time] = process_sino_bwd(self, blobs, n) % Here we extract the individual sinograms from the blobs if (self.using_super_sampling()) @@ -167,14 +153,6 @@ classdef Gt6DVolumeToBlobProjector < Gt6DVolumeProjector sinogram_size, blobs, offs); end bs_time = toc(c); - - if (~isempty(self.shape_functions)) - c = tic(); - sino = self.apply_shape_functions_bp(sino, self.shape_functions{n}); - sf_time = toc(c); - else - sf_time = []; - end end end @@ -236,88 +214,5 @@ classdef Gt6DVolumeToBlobProjector < Gt6DVolumeProjector rethrow(mexc) end end - - function sino_sf = apply_shape_functions_fp(sino, shape_funcs) - sizes_shape_funcs = cat(1, shape_funcs(:).bbsize); - depth_sino_sf = sum(sizes_shape_funcs(:, 3)); - - slice_size = [size(sino, 1), size(sino, 3)]; - sino_sf = zeros(slice_size(1), slice_size(2), depth_sino_sf, 'like', sino); - - sino = permute(sino, [1 3 2]); - -% % Example of same code, but using FFT -% padded_size = slice_size + sizes_shape_funcs(1:2) - 1; -% shifts_sf = [round((padded_size - sizes_shape_funcs(1:2)) / 2), 0]; -% -% padded_sino = zeros([padded_size, size(sino, 3)], 'like', sino); -% shifts_sino = [round((padded_size - slice_size) / 2), 0]; -% padded_sino = gtPlaceSubVolume(padded_sino, sino, shifts_sino); -% padded_sino = fft2(padded_sino); -% -% cumul_pos = 0; -% for ii_b = 1:numel(shape_funcs) -% bl = shape_funcs(ii_b); -% padded_intm = zeros([padded_size, bl.bbsize(3)], 'like', bl.intm); -% padded_intm = gtPlaceSubVolume(padded_intm, bl.intm, shifts_sf); -% padded_intm = fft2(padded_intm); -% -% temp_mult = padded_sino(:, :, ii_b * ones(shape_funcs(ii_b).bbsize(3), 1)) .* padded_intm; -% temp_mult = ifft2(temp_mult); -% inds_u_temp_mult = shifts_sino(1) + (1:slice_size(1)); -% inds_v_temp_mult = shifts_sino(2) + (1:slice_size(2)); -% temp_mult = temp_mult(inds_u_temp_mult, inds_v_temp_mult, :); -% temp_mult = fftshift(real(temp_mult), 1); -% temp_mult = fftshift(temp_mult, 2); -% -% inds_sino_sf = cumul_pos+(1:shape_funcs(ii_b).bbsize(3)); -% sino_sf(:, inds_sino_sf, :) = permute(temp_mult, [1 3 2]); -% -% cumul_pos = cumul_pos + shape_funcs(ii_b).bbsize(3); -% end - - cumul_pos = 0; - for ii_b = 1:numel(shape_funcs) - sino_slice = sino(:, :, ii_b); - sf_intm = shape_funcs(ii_b).intm; - - for ii_w = 1:shape_funcs(ii_b).bbsize(3) - sf_slice = sf_intm(:, :, ii_w); - sino_sf(:, :, cumul_pos+ii_w) = conv2(sino_slice, sf_slice, 'same'); - end - - cumul_pos = cumul_pos + shape_funcs(ii_b).bbsize(3); - end - - sino_sf = permute(sino_sf, [1 3 2]); - end - - function sino = apply_shape_functions_bp(sino_sf, shape_funcs) - num_blobs = numel(shape_funcs); - slice_size = [size(sino_sf, 1), size(sino_sf, 3)]; - sino = zeros(slice_size(1), slice_size(2), num_blobs, 'like', sino_sf); - - sino_sf = permute(sino_sf, [1 3 2]); - - cumul_pos = 0; - for ii_b = 1:num_blobs - sf_w_size = shape_funcs(ii_b).bbsize(3); - - sf_intm = flip(flip(shape_funcs(ii_b).intm, 1), 2); - sino_sf_stack = sino_sf(:, :, cumul_pos+1:cumul_pos+sf_w_size); - - conv_slices = zeros([slice_size, sf_w_size], 'like', sino_sf); - for ii_w = 1:sf_w_size - sino_sf_slice = sino_sf_stack(:, :, ii_w); - sf_slice = sf_intm(:, :, ii_w); - conv_slices(:, :, ii_w) = conv2(sino_sf_slice, sf_slice, 'same'); - end - sino(:, :, ii_b) = sum(conv_slices, 3); - - cumul_pos = cumul_pos + shape_funcs(ii_b).bbsize(3); - end - - sino = permute(sino, [1 3 2]); - end end end \ No newline at end of file