Gt6DVolumeProjector.m

classdef Gt6DVolumeProjector < handle
    properties
        proj_sizes_uv = [];

        volume_ss = 1; % Volume downscaling option
        rspace_ss = 1; % Real-space oversampling
        detector_ss = []; % Detector-wise oversampling

        volume_geometry = {};
        astra_volume_geometry = {};

        geometries = {};
        astra_projection_geometries = {};
        astra_projector_ids = {};

        data_type = 'single';

        num_threads = 1;
        num_gpus = 1;
        jobs_bunch_size = 16;

        statistics = GtTasksStatistics();
    end

    properties (Constant)
        messageNoGPU = 'This machine cannot be used to forward and back project';
    end

    methods (Access = public)
        function self = Gt6DVolumeProjector(vols_size, proj_sizes_uv, varargin)
            % volume geometry (x, y, z)
            self.volume_geometry = vols_size;
            self.astra_volume_geometry = astra_create_vol_geom(vols_size(2), vols_size(1), vols_size(3));

            self.proj_sizes_uv = proj_sizes_uv;

            self = parse_pv_pairs(self, varargin);

            if (exist('astra_mex_direct_c', 'file') ~= 3)
                error('Gt6DVolumeProjector:bad_astra_installation', ...
                    '"astra_mex_direct_c" is not available! Please update or recompile ASTRA')
            end

            num_det = self.get_number_detectors();
            num_geoms = numel(self.geometries{1});
            for ii_d = 2:num_det
                if (num_geoms ~= numel(self.geometries{ii_d}))
                    error('Gt6DVolumeProjector:wrong_argument', ...
                        'All detectors should have the same number of geometries!')
                end
            end

            self.initProjectionGeometry();

            try
                self.num_threads = feature('NumCores');
                self.jobs_bunch_size = self.num_threads;
            catch
            end

            try
                xml_conf = gtConfLoadXML();
                astra_gpus = gtConfGetField(xml_conf, 'astra.gpu');
                for ii_c = 1:numel(astra_gpus)
                    gpus = gtConfFilterAttribute(astra_gpus(ii_c), 'hostname');
                    if (~isempty(gpus))
                        self.num_gpus = gtConfGetField(gpus, 'count');
                        try
                            gpus_indx = gtConfGetField(gpus, 'indx');
                        catch
                            gpus_indx = 0:self.num_gpus-1;
                        end
                        astra_mex('set_gpu_index', gpus_indx);
                        break;
                    end
                end
            catch mexc
                self.num_gpus = 1;
                gtPrintException(mexc, ...
                    'No Astra gpu information, defaulting to 1 gpu only');
            end
        end

        function delete(self)
            self.reset_geometry();
        end

        function initProjectionGeometry(self)
            % Let's clean up previous geometries
            self.reset_geometry();

            % Volume Downscaling option and similar
            opts = struct( ...
                'VoxelSuperSampling', self.volume_ss * self.rspace_ss, ...
                'DetectorSuperSampling', self.rspace_ss, ...
                'GPUindex', -1 );

            num_geoms = numel(self.geometries{1});
            num_det = self.get_number_detectors();

            for ii_d = 1:num_det
                det_ss_factor = ceil(self.detector_ss(ii_d) - 0.1);
                opts_n = opts;
                opts_n.DetectorSuperSampling = opts_n.DetectorSuperSampling * det_ss_factor;

                self.astra_projection_geometries{ii_d} = cell(num_geoms, 1);
                self.astra_projector_ids{ii_d} = cell(num_geoms, 1);

                for n = 1:num_geoms
                    geom = self.geometries{ii_d}{n};

                    self.astra_projection_geometries{ii_d}{n} = astra_create_proj_geom(...
                        'parallel3d_vec', self.proj_sizes_uv(ii_d, 2), self.proj_sizes_uv(ii_d, 1), geom);

                    self.astra_projector_ids{ii_d}{n} = astra_create_projector('cuda3d', ...
                        self.astra_projection_geometries{ii_d}{n}, self.astra_volume_geometry, opts_n);
                end
            end
        end

        function num_det = get_number_detectors(self)
            num_det = size(self.proj_sizes_uv, 1);
        end

        function num_geometries = get_number_geometries(self)
            num_geometries = numel(self.astra_projector_ids{1});
        end

        function chunk_size = get_jobs_chunk_size(self)
            chunk_size = self.num_gpus * self.jobs_bunch_size;
        end

        function printStats(self)
            self.statistics.printStats()
        end

        function stats = get_statistics(self)
            stats = self.statistics;
        end
    end

    methods (Access = protected)
        function sinogram = fwd_project_volumes_to_sinos(self, volume, det_ind, n)
        % Basic Fwd-Projection function

            sinogram = astra_mex_direct_c('FP3D', [self.astra_projector_ids{det_ind}{n}], volume);
        end

        function volume = bwd_project_sinos_to_volumes(self, sinogram, det_ind, n)
        % Basic Bwd-Projection function

            volume = astra_mex_direct_c('BP3D', [self.astra_projector_ids{det_ind}{n}], sinogram);
        end

        function reset_geometry(self)
            num_det = self.get_number_detectors();
            self.astra_projection_geometries = cell(num_det, 1);

            if (~isempty(self.astra_projector_ids))
                for ii_d = 1:numel(self.astra_projector_ids)
                    if (~isempty(self.astra_projector_ids{ii_d}))
                        astra_mex_projector3d('delete', self.astra_projector_ids{ii_d}{:});
                    end
                end
                self.astra_projector_ids = cell(num_det, 1);
            end
        end
    end
end