classdef GtGrainODFuvwSolver < GtGrainODFwSolver
properties
size_rspace_volume = [];
voxel_centers = [];
rspace_oversampling = 1;
rspace_downscaling = 1;
end
methods (Access = public)
function self = GtGrainODFuvwSolver(parameters, varargin)
self = self@GtGrainODFwSolver(parameters, varargin{:});
end
end
methods (Access = public) % Low Level API
function build_orientation_sampling(self, ref_gr, det_index)
self.build_orientation_sampling@GtGrainODFwSolver(ref_gr, det_index);
self.chose_voxel_centers();
num_vox_centers = size(self.voxel_centers, 1);
recgeo = self.parameters.recgeo(self.sampler.detector_index);
switch (self.shape_functions_type)
case 'nw2uvw'
sf_nw = gtDefShapeFunctionsFwdProj(self.sampler, ...
'shape_function_type', 'nw', ...
'data_type', self.data_type);
% sf_nw = gtDefShapeFunctionsCreateNW(self.sampler, 2, ...
% 'data_type', self.data_type);
for ii_v = 1:num_vox_centers
fprintf('(%03d/%03d) ', ii_v, num_vox_centers);
disp_sampler = self.sampler.regenerate_displaced(self.voxel_centers(ii_v, :));
self.shape_functions{ii_v} = gtDefShapeFunctionsNW2UVW( ...
disp_sampler, sf_nw, 'recenter_sf', false, ...
'rspace_oversampling', self.rspace_downscaling, ...
'rspace_voxel_size', recgeo.voxsize * self.rspace_downscaling);
end
case 'uvw'
shape_functions_name = sprintf('tmp_shape_functions_uvw_%d_%d', ...
size(self.voxel_centers, 1), self.sampler.get_total_num_orientations());
if (evalin('base', ['exist(''' shape_functions_name ''', ''var'')']))
fprintf('Loading shape functions..')
c = tic();
self.shape_functions = evalin('base', shape_functions_name);
fprintf('\b\b: Done in %g seconds\n', toc(c));
elseif (~exist([shape_functions_name '.mat'], 'file'))
for ii_v = 1:num_vox_centers
fprintf('(%3d/%03d) ', ii_v, num_vox_centers);
disp_sampler = self.sampler.regenerate_displaced(self.voxel_centers(ii_v, :));
self.shape_functions{ii_v} = gtDefShapeFunctionsFwdProj( ...
disp_sampler, 'recenter_sf', false, ...
'data_type', self.data_type, ...
'rspace_voxel_size', recgeo.voxsize * self.rspace_downscaling);
end
sfs = self.shape_functions;
save([shape_functions_name '.mat'], 'sfs', '-v7.3')
assignin('base', shape_functions_name, sfs);
else
fprintf('Loading shape functions..')
c = tic();
sfs_file = load([shape_functions_name '.mat']);
self.shape_functions = sfs_file.sfs;
assignin('base', shape_functions_name, sfs_file.sfs);
fprintf('\b\b: Done in %g seconds\n', toc(c));
end
end
end
function build_sinogram(self)
if (self.verbose)
fprintf('Building sinogram..')
c = tic();
end
ref_sel = self.sampler.selected;
bls = self.sampler.ref_gr.proj(self.sampler.detector_index).bl(ref_sel);
ref_inc = self.sampler.ondet(self.sampler.included);
ref_gr = self.sampler.get_reference_grain();
ref_ws = ref_gr.allblobs(self.sampler.detector_index).detector.uvw(ref_inc(ref_sel), 3);
num_blobs = numel(bls);
blob_sizes = cat(1, bls(:).bbsize);
blob_u_lims = cat(1, bls(:).bbuim);
blob_v_lims = cat(1, bls(:).bbvim);
blob_w_lims = cat(1, bls(:).bbwim);
blob_uvw_lims = permute(cat(3, blob_u_lims, blob_v_lims, blob_w_lims), [1 3 2]);
with_shape_functions = ~isempty(self.shape_functions);
% Building W information
if (with_shape_functions)
if (any(strcmpi(self.shape_functions_type, {'uvw', 'nw2uvw'})))
for ii_o = numel(self.shape_functions{1}):-1:1
proj_w_lims(:, :, ii_o) = cat(1, self.shape_functions{1}{ii_o}(:).bbwim);
end
else
for ii_o = numel(self.shape_functions):-1:1
proj_w_lims(:, :, ii_o) = cat(1, self.shape_functions{ii_o}(:).bbwim);
end
end
proj_w_mins = reshape(proj_w_lims(:, 1, :), num_blobs, []);
proj_w_maxs = reshape(proj_w_lims(:, 2, :), num_blobs, []);
proj_w_mins = gtGrainAnglesTabularFix360deg(proj_w_mins, ref_ws, self.parameters);
proj_w_maxs = gtGrainAnglesTabularFix360deg(proj_w_maxs, ref_ws, self.parameters);
proj_w_lims = [min(proj_w_mins, [], 2), max(proj_w_maxs, [], 2)];
delta_ws = proj_w_lims(:, 2) - proj_w_lims(:, 1) + 1;
else
[delta_ws, proj_w_lims] = self.sampler.get_omega_deviations(with_shape_functions);
delta_ws = delta_ws(self.sampler.selected);
proj_w_lims = proj_w_lims(self.sampler.selected, :);
end
num_rspace_voxels = size(self.voxel_centers, 1);
% Building UV information
if (with_shape_functions && any(strcmpi(self.shape_functions_type, {'uvw', 'nw2uvw'})))
for ii_v = num_rspace_voxels:-1:1
for ii_o = numel(self.shape_functions{1}):-1:1
proj_u_lims(:, :, ii_o, ii_v) = cat(1, self.shape_functions{ii_v}{ii_o}(:).bbuim);
proj_v_lims(:, :, ii_o, ii_v) = cat(1, self.shape_functions{ii_v}{ii_o}(:).bbvim);
end
end
% By expressing the third index with the column :, and
% not expressing the fourth at all, we do an implicit
% reshape from 4 imentions to 3
proj_u_lims = [min(proj_u_lims(:, 1, :), [], 3), max(proj_u_lims(:, 2, :), [], 3)];
proj_v_lims = [min(proj_v_lims(:, 1, :), [], 3), max(proj_v_lims(:, 2, :), [], 3)];
else
for ii_v = num_rspace_voxels:-1:1
disp_sampler = self.sampler.regenerate_displaced(self.voxel_centers(ii_v, :));
[~, proj_uv_lims(:, :, :, ii_v)] = disp_sampler.get_uv_deviations();
end
proj_uv_lims = proj_uv_lims(self.sampler.selected, :, :);
proj_u_lims = [min(proj_uv_lims(:, 1, 1, :), [], 4), max(proj_uv_lims(:, 2, 1, :), [], 4)];
proj_v_lims = [min(proj_uv_lims(:, 1, 2, :), [], 4), max(proj_uv_lims(:, 2, 2, :), [], 4)];
end
proj_uv_lims = cat(3, proj_u_lims, proj_v_lims);
delta_uvs = reshape(proj_uv_lims(:, 2, :) - proj_uv_lims(:, 1, :), [], 2);
chosen_sizes = max(blob_sizes, [delta_uvs, delta_ws]);
proj_uvw_lims = permute(cat(3, proj_uv_lims, proj_w_lims), [1 3 2]);
num_uvws = max(chosen_sizes, [], 1) + 2;
self.size_sino = [num_uvws, num_blobs];
self.sino = zeros(self.size_sino, self.data_type);
% In the case that the sampled orientations determine the
% shift, we have to take it into account!
align_uvw_shift = blob_uvw_lims(:, :, 1) - proj_uvw_lims(:, :, 1);
align_uvw_shift(align_uvw_shift < 0) = 0;
self.pre_paddings = floor((num_uvws(ones(num_blobs, 1), :) - chosen_sizes) / 2) + 1 + align_uvw_shift;
for ii_b = 1:num_blobs
ints_u_interval = self.pre_paddings(ii_b, 1):(self.pre_paddings(ii_b, 1) + blob_sizes(ii_b, 1) -1);
ints_v_interval = self.pre_paddings(ii_b, 2):(self.pre_paddings(ii_b, 2) + blob_sizes(ii_b, 2) -1);
ints_w_interval = self.pre_paddings(ii_b, 3):(self.pre_paddings(ii_b, 3) + blob_sizes(ii_b, 3) -1);
masked_blob = bls(ii_b).intm;
masked_blob(~bls(ii_b).mask) = 0;
self.sino(ints_u_interval, ints_v_interval, ints_w_interval, ii_b) = masked_blob;
end
self.sino = reshape(self.sino, [], 1);
if (self.verbose)
fprintf('\b\b: Done in %g seconds.\n - Sino size: [%d, %d, %d, %d]\n', toc(c), self.size_sino);
end
end
function chose_voxel_centers(self)
recgeo = self.parameters.recgeo(self.sampler.detector_index);
ref_gr = self.sampler.get_reference_grain();
proj = ref_gr.proj(self.sampler.detector_index);
rspace_vol_dims = ceil([proj.vol_size_y, proj.vol_size_x, proj.vol_size_z] / self.rspace_downscaling);
max_dists_from_center = (rspace_vol_dims - 1 ./ self.rspace_oversampling) / 2;
self.size_rspace_volume = rspace_vol_dims .* self.rspace_oversampling;
tot_voxels = prod(self.size_rspace_volume);
self.voxel_centers = zeros(tot_voxels, 3);
counter = 1;
for ss_z = linspace(-max_dists_from_center(3), max_dists_from_center(3), self.size_rspace_volume(3))
for ss_y = linspace(-max_dists_from_center(2), max_dists_from_center(2), self.size_rspace_volume(2))
for ss_x = linspace(-max_dists_from_center(1), max_dists_from_center(1), self.size_rspace_volume(1))
offset_voxel = [ss_x, ss_y, ss_z] .* recgeo.voxsize * self.rspace_downscaling;
self.voxel_centers(counter, :) = ref_gr.center + offset_voxel;
counter = counter + 1;
end
end
end
end
function build_projection_matrices_sf_none(self)
ref_sel = self.sampler.selected;
det_ind = self.sampler.detector_index;
detgeo = self.parameters.detgeo(det_ind);
labgeo = self.parameters.labgeo;
samgeo = self.parameters.samgeo;
om_step = gtAcqGetOmegaStep(self.parameters, det_ind);
bls = self.sampler.ref_gr.proj(det_ind).bl(ref_sel);
num_uvws = self.get_num_uvws();
num_blobs = self.get_num_blobs();
bls_bbus = cat(1, bls(:).bbuim);
bls_bbvs = cat(1, bls(:).bbvim);
bls_bbws = cat(1, bls(:).bbwim);
min_uvw_conds = [bls_bbus(:, 1), bls_bbvs(:, 1), bls_bbws(:, 1)] - self.pre_paddings + 1;
max_uvw_conds = min_uvw_conds + num_uvws(ones(num_blobs, 1), :) - 1;
grid_gr = self.sampler.get_orientations();
num_orients = numel(grid_gr);
num_voxels = size(self.voxel_centers, 1);
self.S = cell(num_voxels, 1);
self.St = cell(num_voxels, 1);
rot_l2s = cell(num_orients, 1);
for ii_o = 1:num_orients
rot_l2s{ii_o} = grid_gr{ii_o}.allblobs(det_ind).srot(:, :, ref_sel);
end
rot_l2s = cat(3, rot_l2s{:});
rot_s2l = permute(rot_l2s, [2 1 3]);
dveclab = cell(num_orients, 1);
for ii_o = 1:num_orients
dveclab{ii_o} = grid_gr{ii_o}.allblobs(det_ind).dvec(ref_sel, :);
end
dveclab = cat(1, dveclab{:});
for ii_v = 1:num_voxels
b_us = [];
b_vs = [];
b_ws = [];
b_cs = [];
b_is = [];
b_os = [];
num_chars = fprintf('%03d/%03d', ii_v, num_voxels);
gcsam_v = self.voxel_centers(ii_v, :);
gcsam_v = gcsam_v(ones(1, size(rot_s2l, 3)), :);
gclab_v = gtGeoSam2Lab(gcsam_v, rot_s2l, labgeo, samgeo, false);
uv_tot = gtFedPredictUVMultiple([], dveclab', gclab_v', ...
detgeo.detrefpos', detgeo.detnorm', detgeo.Qdet, ...
[detgeo.detrefu, detgeo.detrefv]');
uv_tot = reshape(uv_tot, 2, [], num_orients);
for ii_o = 1:num_orients
% W part
ws = grid_gr{ii_o}.allblobs(det_ind).omega(ref_sel) / om_step;
% UV part
uvw = [uv_tot(:, :, ii_o)', ws];
num_pos = size(uvw, 1);
min_uvws = floor(uvw);
max_uvws = min_uvws + 1;
max_uvw_cs = uvw - min_uvws;
min_uvw_cs = 1 - max_uvw_cs;
% Acceptance criteria
ok_uvw_mins = (min_uvws >= min_uvw_conds) & (min_uvws <= max_uvw_conds);
ok_uvw_maxs = (max_uvws >= min_uvw_conds) & (max_uvws <= max_uvw_conds) & (max_uvw_cs > eps('single'));
u_oks = [ok_uvw_mins(:, 1), ok_uvw_maxs(:, 1)];
v_oks = [ok_uvw_mins(:, 2), ok_uvw_maxs(:, 2)];
w_oks = [ok_uvw_mins(:, 3), ok_uvw_maxs(:, 3)];
u_cs = [min_uvw_cs(:, 1), max_uvw_cs(:, 1)];
v_cs = [min_uvw_cs(:, 2), max_uvw_cs(:, 2)];
w_cs = [min_uvw_cs(:, 3), max_uvw_cs(:, 3)];
pos_us = [min_uvws(:, 1), max_uvws(:, 1)] - bls_bbus(:, [1 1]) + self.pre_paddings(:, [1 1]);
pos_vs = [min_uvws(:, 2), max_uvws(:, 2)] - bls_bbvs(:, [1 1]) + self.pre_paddings(:, [2 2]);
pos_ws = [min_uvws(:, 3), max_uvws(:, 3)] - bls_bbws(:, [1 1]) + self.pre_paddings(:, [3 3]);
uvw_oks = u_oks(:, [1 1 1 1 2 2 2 2]) & v_oks(:, [1 1 2 2 1 1 2 2]) & w_oks(:, [1 2 1 2 1 2 1 2]);
uvw_cs = u_cs(:, [1 1 1 1 2 2 2 2]) .* v_cs(:, [1 1 2 2 1 1 2 2]) .* w_cs(:, [1 2 1 2 1 2 1 2]);
pos_us = pos_us(:, [1 1 1 1 2 2 2 2]);
pos_vs = pos_vs(:, [1 1 2 2 1 1 2 2]);
pos_ws = pos_ws(:, [1 2 1 2 1 2 1 2]);
indx = find(uvw_oks);
indx_mod = mod(indx - 1, num_pos) + 1;
b_us = [b_us; pos_us(indx)];
b_vs = [b_vs; pos_vs(indx)];
b_ws = [b_ws; pos_ws(indx)];
b_cs = [b_cs; uvw_cs(indx)];
b_is = [b_is; indx_mod];
b_os = [b_os; ii_o(ones(numel(indx), 1))];
end
% sino_indx = sub2ind(self.size_sino, b_us, b_vs, b_ws, b_is);
sino_indx = b_us ...
+ num_uvws(1) * (b_vs - 1 ...
+ num_uvws(2) * (b_ws - 1 ...
+ num_uvws(3) * (b_is - 1)));
self.S{ii_v} = sparse( ...
sino_indx, b_os, b_cs, ...
numel(self.sino), num_orients);
self.St{ii_v} = self.S{ii_v}';
fprintf(repmat('\b', [1 num_chars]));
end
end
function build_projection_matrices_sf_w(self)
ref_sel = self.sampler.selected;
det_ind = self.sampler.detector_index;
detgeo = self.parameters.detgeo(det_ind);
labgeo = self.parameters.labgeo;
samgeo = self.parameters.samgeo;
bls = self.sampler.ref_gr.proj(self.sampler.detector_index).bl(ref_sel);
num_uvws = self.get_num_uvws();
num_blobs = self.get_num_blobs();
bls_bbws = cat(1, bls(:).bbwim);
bls_bbus = cat(1, bls(:).bbuim);
bls_bbvs = cat(1, bls(:).bbvim);
min_uvw_conds = [bls_bbus(:, 1), bls_bbvs(:, 1), bls_bbws(:, 1)] - self.pre_paddings + 1;
max_uvw_conds = min_uvw_conds + num_uvws(ones(num_blobs, 1), :) - 1;
grid_gr = self.sampler.get_orientations();
num_orients = numel(grid_gr);
num_voxels = size(self.voxel_centers, 1);
self.S = cell(num_voxels, 1);
self.St = cell(num_voxels, 1);
% Extracting orientations' W shape functions
oris_w_cs = cell(num_orients, 1);
oris_ws = cell(num_orients, 1);
oris_w_is = cell(num_orients, 1);
for ii_o = 1:num_orients
sf = self.shape_functions{ii_o};
w_cs = cat(1, sf(:).intm);
oris_w_cs{ii_o} = reshape(w_cs, [], 1);
ws_sizes = [sf(:).bbsize];
bbwims = cat(1, sf(:).bbwim);
ws = cell(num_blobs, 1);
w_is = cell(num_blobs, 1);
for ii_b = 1:num_blobs
ws{ii_b} = bbwims(ii_b, 1):bbwims(ii_b, 2);
w_is{ii_b} = ii_b(ones(ws_sizes(ii_b), 1));
end
oris_ws{ii_o} = reshape([ws{:}], [], 1);
oris_w_is{ii_o} = cat(1, w_is{:});
end
% Preparing for computing the UV deviations due to the XYZ
% positions of the real-space voxels
rot_l2s = cell(num_orients, 1);
for ii_o = 1:num_orients
rot_l2s{ii_o} = grid_gr{ii_o}.allblobs(det_ind).srot(:, :, ref_sel);
end
rot_l2s = cat(3, rot_l2s{:});
rot_s2l = permute(rot_l2s, [2 1 3]);
dveclab = cell(num_orients, 1);
for ii_o = 1:num_orients
dveclab{ii_o} = grid_gr{ii_o}.allblobs(det_ind).dvec(ref_sel, :);
end
dveclab = cat(1, dveclab{:});
for ii_v = 1:num_voxels
b_us = [];
b_vs = [];
b_ws = [];
b_cs = [];
b_is = [];
b_os = [];
num_chars = fprintf('%03d/%03d', ii_v, num_voxels);
gcsam_v = self.voxel_centers(ii_v, :);
gcsam_v = gcsam_v(ones(1, size(rot_s2l, 3)), :);
gclab_v = gtGeoSam2Lab(gcsam_v, rot_s2l, labgeo, samgeo, false);
uv_tot = gtFedPredictUVMultiple([], dveclab', gclab_v', ...
detgeo.detrefpos', detgeo.detnorm', detgeo.Qdet, ...
[detgeo.detrefu, detgeo.detrefv]');
uv_tot = reshape(uv_tot, 2, [], num_orients);
for ii_o = 1:num_orients
w_cs = oris_w_cs{ii_o};
ws = oris_ws{ii_o};
w_is = oris_w_is{ii_o};
% UV part
uv = uv_tot(:, :, ii_o)';
min_uvs = floor(uv);
max_uvs = min_uvs + 1;
max_uv_cs = uv - min_uvs;
min_uv_cs = 1 - max_uv_cs;
% Acceptance criteria
ok_uv_mins = (min_uvs >= min_uvw_conds(:, 1:2)) & (min_uvs <= max_uvw_conds(:, 1:2)) & (min_uv_cs > eps('single'));
ok_uv_maxs = (max_uvs >= min_uvw_conds(:, 1:2)) & (max_uvs <= max_uvw_conds(:, 1:2)) & (max_uv_cs > eps('single'));
u_oks = [ok_uv_mins(:, 1), ok_uv_maxs(:, 1)];
v_oks = [ok_uv_mins(:, 2), ok_uv_maxs(:, 2)];
u_cs = [min_uv_cs(:, 1), max_uv_cs(:, 1)];
v_cs = [min_uv_cs(:, 2), max_uv_cs(:, 2)];
pos_us = [(min_uvs(:, 1) - min_uvw_conds(:, 1) + 1), (max_uvs(:, 1) - min_uvw_conds(:, 1) + 1)];
pos_vs = [(min_uvs(:, 2) - min_uvw_conds(:, 2) + 1), (max_uvs(:, 2) - min_uvw_conds(:, 2) + 1)];
pos_ws = ws - min_uvw_conds(w_is, 3) + 1;
uvw_oks = u_oks(w_is, [1 1 2 2]) & v_oks(w_is, [1 2 1 2]);
uvw_cs = u_cs(w_is, [1 1 2 2]) .* v_cs(w_is, [1 2 1 2]) .* w_cs(:, [1 1 1 1]);
pos_us = pos_us(w_is, [1 1 2 2]);
pos_vs = pos_vs(w_is, [1 2 1 2]);
pos_ws = pos_ws(:, [1 1 1 1]);
indx = find(uvw_oks);
indx_mod = mod(indx - 1, numel(w_is)) + 1;
wrong_w = ws < min_uvw_conds(w_is, 3) | ws > max_uvw_conds(w_is, 3);
if (any(wrong_w))
ii_o
find(wrong_w)
[min_uvw_conds(w_is(wrong_w), 3), ws(wrong_w), max_uvw_conds(w_is(wrong_w), 3), ...
(max_uvw_conds(w_is(wrong_w), 3) - min_uvw_conds(w_is(wrong_w), 3) + 1), ...
(ws(wrong_w) - min_uvw_conds(w_is(wrong_w), 3) + 1) ]
end
b_us = [b_us; pos_us(indx)];
b_vs = [b_vs; pos_vs(indx)];
b_ws = [b_ws; pos_ws(indx)];
b_cs = [b_cs; uvw_cs(indx)];
b_is = [b_is; w_is(indx_mod)];
b_os = [b_os; ii_o(ones(numel(indx), 1))];
end
% sino_indx = sub2ind(self.size_sino, b_us, b_vs, b_ws, b_is);
sino_indx = b_us ...
+ self.size_sino(1) * (b_vs - 1 ...
+ self.size_sino(2) * (b_ws - 1 ...
+ self.size_sino(3) * (b_is - 1)));
self.S{ii_v} = sparse( ...
sino_indx, b_os, b_cs, ...
numel(self.sino), num_orients);
self.St{ii_v} = self.S{ii_v}';
fprintf(repmat('\b', [1 num_chars]));
end
end
function build_projection_matrices_sf_uvw(self)
ref_sel = self.sampler.selected;
bls = self.sampler.ref_gr.proj(self.sampler.detector_index).bl(ref_sel);
num_blobs = self.get_num_blobs();
num_ws = self.get_num_ws();
bls_bbws = cat(1, bls(:).bbwim);
min_w_conds = bls_bbws(:, 1) - self.pre_paddings(:, 3) + 1;
max_w_conds = min_w_conds + num_ws - 1;
bls_bbus = cat(1, bls(:).bbuim);
bls_bbvs = cat(1, bls(:).bbvim);
grid_gr = self.sampler.get_orientations();
num_orients = numel(grid_gr);
num_voxels = size(self.voxel_centers, 1);
self.S = cell(num_voxels, 1);
self.St = cell(num_voxels, 1);
bls_uvw_orig = cat(1, bls(:).bbuim, bls(:).bbvim, bls(:).bbwim);
bls_uvw_orig = reshape(bls_uvw_orig, [], 3, 2);
bls_uvw_orig = bls_uvw_orig(:, :, 1) - self.pre_paddings + 1;
for ii_v = 1:num_voxels
b_us = [];
b_vs = [];
b_ws = [];
b_cs = [];
b_is = [];
b_os = [];
num_chars = fprintf('%03d/%03d', ii_v, num_voxels);
for ii_o = 1:num_orients
sf = self.shape_functions{ii_v}{ii_o};
sf_uvw_shift = cat(1, sf(:).bbuim, sf(:).bbvim, sf(:).bbwim);
sf_uvw_shift = reshape(sf_uvw_shift, [], 3, 2);
sf_uvw_shift = sf_uvw_shift(:, :, 1) - bls_uvw_orig;
uvw_cs = cell(num_blobs, 1);
bl_is = cell(num_blobs, 1);
bl_shifts = cell(num_blobs, 1);
for ii_b = 1:num_blobs
inds = find(sf(ii_b).intm);
uvw_cs{ii_b} = sf(ii_b).intm(inds);
% [inds_u, inds_v, inds_w] = ind2sub(sf(ii_b).bbsize, inds);
inds_u = mod(inds, sf(ii_b).bbsize(1));
inds_v = mod(floor(inds / sf(ii_b).bbsize(1)) + 1, sf(ii_b).bbsize(2));
inds_w = floor(inds / prod(sf(ii_b).bbsize(1:2))) + 1;
bl_is{ii_b} = ii_b(ones(numel(inds), 1));
bl_shifts{ii_b} = [inds_u, inds_v, inds_w] + sf_uvw_shift(bl_is{ii_b}, :);
end
uvw_cs = cat(1, uvw_cs{:});
bl_is = cat(1, bl_is{:});
bl_shifts = cat(1, bl_shifts{:});
b_us = [b_us; bl_shifts(:, 1)];
b_vs = [b_vs; bl_shifts(:, 2)];
b_ws = [b_ws; bl_shifts(:, 3)];
b_cs = [b_cs; uvw_cs];
b_is = [b_is; bl_is];
b_os = [b_os; ii_o(ones(numel(bl_is), 1))];
end
% sino_indx = sub2ind(self.size_sino, b_us, b_vs, b_ws, b_is);
sino_indx = b_us ...
+ self.size_sino(1) * (b_vs - 1 ...
+ self.size_sino(2) * (b_ws - 1 ...
+ self.size_sino(3) * (b_is - 1)));
self.S{ii_v} = sparse( ...
sino_indx, b_os, double(b_cs), ...
numel(self.sino), num_orients);
self.St{ii_v} = self.S{ii_v}';
fprintf(repmat('\b', [1 num_chars]));
end
end
function vol = get_volume(self)
vol = reshape(sum(self.volume, 2), self.size_volume);
end
function vol = get_volume_6D(self)
vol = reshape(self.volume, [self.size_volume, self.size_rspace_volume]);
end
end
methods (Access = protected)
function y = fp(self, x)
y = self.S{1} * x(:, 1);
for ii_s = 2:numel(self.S)
y = y + self.S{ii_s} * x(:, ii_s);
end
end
function y = bp(self, x)
for ii_s = numel(self.St):-1:1
y(:, ii_s) = self.St{ii_s} * x;
end
end
function num_uvws = get_num_uvws(self)
num_uvws = self.size_sino(1:3);
end
function num_ws = get_num_ws(self)
num_ws = self.size_sino(3);
end
function num_bs = get_num_blobs(self)
num_bs = self.size_sino(4);
end
end
end