GtGrainODFuvwSolver.m

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