Gt6DVolumeToBlobProjector.m

classdef Gt6DVolumeToBlobProjector < Gt6DVolumeProjector
    properties
        % Data
        blobs;

        offsets = {};

        sinogram_sizes = zeros(0, 3, 0);

        % Weights
        fwd_weights = {};
        bwd_weights = {};

        blob_sinogram_c_functions = true;
    end

    methods (Access = public)
        function self = Gt6DVolumeToBlobProjector(vols_size, blobs, proj_sizes_uv, varargin)
            detector_ss = cat(1, blobs(:).det_ss);
            args = [varargin, {'detector_ss', detector_ss}];

            self = self@Gt6DVolumeProjector(vols_size, proj_sizes_uv, args{:});

            self.blobs = blobs;

            have_multi_sino_functions = ...
                (exist('gtCxx6DMultipleSinosToBlobs', 'file') == 3) ...
                && (exist('gtCxx6DBlobsToMultipleSinos', 'file') == 3);
            if (~have_multi_sino_functions)
                warning('Gt6DVolumeToBlobProjector:missing_fast_functions', ...
                    '"gtCxx6DMultipleSinosToBlobs" or/and "gtCxx6DBlobsToMultipleSinos" not available! Please recompile mex files')
                self.blob_sinogram_c_functions = false;
            end

            num_orients = self.get_number_geometries();
            ones_n_ors = ones(num_orients, 1);

            num_det = self.get_number_detectors();
            self.sinogram_sizes = zeros(num_orients, 3, num_det);
            for ii_d = 1:num_det
                geoms = self.geometries{ii_d};

                self.sinogram_sizes(:, :, ii_d) = [ ...
                    self.proj_sizes_uv(ii_d * ones_n_ors, 1), zeros(num_orients, 1), self.proj_sizes_uv(ii_d * ones_n_ors, 2)];
                for ii_o = 1:num_orients
                    self.sinogram_sizes(ii_o, 2, ii_d) = size(geoms{ii_o}, 1);
                end
            end
        end
    end

    methods (Access = protected)
        function [blobs, fp_time, bs_time, pbs_time] = fwd_project_volume(self, blobs, volume, det_ind, inds)
        % 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
        % output and recompose the slices into blobs

            c = tic();
            sino = self.fwd_project_volumes_to_sinos(volume, det_ind, inds);
            fp_time = toc(c);

            c = tic();
            [blobs, pbs_time] = self.process_sinos_fwd(blobs, sino, det_ind, inds);
            bs_time = toc(c);
        end

        function [volumes, bp_time, bs_time, pbs_time] = bwd_project_volume(self, blobs, det_ind, inds)
        % 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

            c = tic();
            [sinos, pbs_time] = self.process_sinos_bwd(blobs, det_ind, inds);
            bs_time = toc(c);

            c = tic();
            volumes = self.bwd_project_sinos_to_volumes(sinos, det_ind, inds);
            bp_time = toc(c);
        end

        function [blobs, bs_time] = process_sinos_fwd(self, blobs, sinos, det_ind, ns)
        % Here we add the individual sinograms to the blobs

            offs = self.offsets{det_ind}(ns);

            c = tic();

            if (self.algo_ops_c_functions)
                try
                    blobs = gtCxx6DMultipleSinosToBlobs(blobs, sinos, offs, 'threads', self.num_threads);
                catch mexc
                    fprintf('Error raised by one of orientations (detector: %d):%s\n', det_ind, sprintf(' %d', ns))
                    rethrow(mexc)
                end
            else
                try
                    for ii_o = 1:numel(inds)
                        blobs = self.single_sinogram_to_blobs(blobs, sinos{ii_o}, offs{ii_o});
                    end
                catch mexc
                    fprintf('Error raised by orientation (detector: %d): %d\n', det_ind, n)
                    rethrow(mexc)
                end
            end
            bs_time = toc(c);
        end

        function [sinos, bs_time] = process_sinos_bwd(self, blobs, det_ind, ns)
        % Here we extract the individual sinograms from the blobs

            offs = self.offsets{det_ind}(ns);

            sinos_size = self.sinogram_sizes(ns, :, det_ind);

            c = tic();
            if (self.algo_ops_c_functions)
                try
                    sinos = gtCxx6DBlobsToMultipleSinos(sinos_size, blobs, offs, 'threads', self.num_threads);
                catch mexc
                    fprintf('Error raised by one of orientations (detector: %d):%s\n', det_ind, sprintf(' %d', ns))
                    rethrow(mexc)
                end
            else
                sinos = cell(num_orients, 1);
                for ii_o = 1:num_orients
                    try
                        sinos{ii_o} = self.blobs_to_single_sinogram(sinos_size(ii_o, :), blobs, offs{det_ind}{ii_o});
                    catch mexc
                        fprintf('Error raised by orientation (detector: %d): %d\n', det_ind, ns)
                        rethrow(mexc)
                    end
                end
            end
            bs_time = toc(c);
        end
    end

    methods (Access = public, Static)
        function blobs = single_sinogram_to_blobs(blobs, sino, coeffs)
            try
                sino_offsets = coeffs.sino_offsets(:);
                blob_offsets = coeffs.blob_offsets(:);
                proj_offsets = coeffs.proj_offsets(:);
                proj_coeffs = coeffs.proj_coeffs(:);

                if (isfield(coeffs, 'proj_shifts_uv'))
                    shifts_uv = coeffs.proj_shifts_uv;
                else
                    shifts_uv = zeros(2, numel(proj_offsets));
                end

                for m = 1:numel(proj_offsets)
                    p = proj_offsets(m);
                    b = blob_offsets(m);
                    c = proj_coeffs(m);
                    s = sino_offsets(m);

                    blobs_us = (1:size(sino, 1)) + shifts_uv(1, m);
                    blobs_vs = (1:size(sino, 3)) + shifts_uv(2, m);

                    blobs{b}(blobs_us, p, blobs_vs) = blobs{b}(blobs_us, p, blobs_vs) + c * sino(:, s, :);
%                     blobs{b}(:, p, :) = blobs{b}(:, p, :) + c * sino(:, s, :);
                end
            catch mexc
                fprintf('\n\nError!!\n')
                fprintf('  indices: b %d, p %d, s %d, m %d\n', ...
                    b, p, s, m)
                fprintf('  size blob (%d %d %d), size sino (%d %d %d)\n\n', ...
                    size(blobs{b}), size(sino))
                fprintf('  Limits over p: [%d, %d]\n', ...
                    min(proj_offsets), max(proj_offsets))
                fprintf('  Limits over s: [%d, %d]\n\n', ...
                    min(sino_offsets), max(sino_offsets))
                rethrow(mexc)
            end
        end

        function sino = blobs_to_single_sinogram(sinogram_size, blobs, coeffs)
            try
                sino = zeros(sinogram_size, 'single');

                sino_offsets = coeffs.sino_offsets(:);
                blob_offsets = coeffs.blob_offsets(:);
                proj_offsets = coeffs.proj_offsets(:);
                proj_coeffs = coeffs.proj_coeffs(:);

                if (isfield(coeffs, 'proj_shifts_uv'))
                    shifts_uv = coeffs.proj_shifts_uv;
                else
                    shifts_uv = zeros(2, numel(proj_offsets));
                end

                for m = 1:numel(proj_offsets)
                    p = proj_offsets(m);
                    b = blob_offsets(m);
                    c = proj_coeffs(m);
                    s = sino_offsets(m);

                    blobs_us = (1:sinogram_size(1)) + shifts_uv(1, m);
                    blobs_vs = (1:sinogram_size(3)) + shifts_uv(2, m);

                    sino(:, s, :) = sino(:, s, :) + c * blobs{b}(blobs_us, p, blobs_vs);
%                     sino(:, s, :) = sino(:, s, :) + c * blobs{b}(:, p, :);
                end
            catch mexc
                fprintf('\n\nError!!\n')
                fprintf('  indices: b %d, p %d, s %d, m %d\n', ...
                    b, p, s, m)
                fprintf('  size blob (%d %d %d), size sino (%d %d %d)\n\n', ...
                    size(blobs{b}), size(sino))
                fprintf('  Limits over p: [%d, %d]\n', ...
                    min(proj_offsets), max(proj_offsets))
                fprintf('  Limits over s: [%d, %d]\n\n', ...
                    min(sino_offsets), max(sino_offsets))
                rethrow(mexc)
            end
        end
    end
end