Gt6DReconstructionAlgorithmFactory.m

classdef Gt6DReconstructionAlgorithmFactory < handle
    properties
        data_type = 'single';

        det_index = 1;

        volume_downscaling = 1;
        rspace_super_sampling = 1;
        rspace_oversize = 1.2;

        use_predicted_scatter_ints = false;
        use_matrix_row_rescaling = false;

        % Different types of shape functions:
        %   'none'   : no shape functions are used
        %   'w'      : Only omega is considered
        %   'nw'     : (eta, omega) shape functions
        shape_functions_type = 'none';

        % Different types of num_interp:
        %   'flat'    : all blobs are inteprolated with the num_interp
        %   'lorentz' : num_interp is rescaled by the Lorentz factor
        num_interp_type = 'lorentz';

        parameters;
    end

    methods (Access = public)
        function self = Gt6DReconstructionAlgorithmFactory(parameters, varargin)
            self = self@handle();

            self = parse_pv_pairs(self, varargin);

            self.parameters = parameters;
        end

        function [algo, blobs] = get6DAlgorithmSingleRegion(self, sampler, num_interp, varargin)
            % Build Empty volumes
            ref_gr = sampler.get_reference_grain();
            volume_size = self.get_volume_size(ref_gr, self.det_index(1));
            vol_size_mm = volume_size .* self.parameters.recgeo(self.det_index(1)).voxsize * self.volume_downscaling;

            or_groups = [1, sampler.get_total_num_orientations()];

            num_det = numel(self.det_index);

            geometries = cell(num_det, 1);
            offsets = cell(num_det, 1);

            blobs = gt6DRecBlobsDefinition(num_det);
            proj_uv_size = zeros(num_det, 2);

            lambda_det = ones(num_det, 1);
%             lambda_det(2:end) = 2e-1;

            for ii_d = 1:num_det
                det_ind = self.det_index(ii_d);
                fprintf('\nProcessing detector: %d\n', det_ind)

                % Deciding the type of num_interp type
                num_interp_d = self.compute_num_interp(num_interp, sampler, det_ind);

                blobs_w_interp = self.compute_blobs_w_interp(sampler, num_interp_d, det_ind);

                shape_funcs = self.get_shape_functions(sampler, num_interp_d, blobs_w_interp, det_ind);

                algo_params = self.get_algorithm_params(sampler, blobs_w_interp, shape_funcs, vol_size_mm, det_ind);

                blobs(ii_d) = algo_params.blobs;
                if (isempty(blobs(ii_d).psf) ...
                        && isfield(self.parameters.rec.grains.options, 'psf') ...
                        && ~isempty(self.parameters.rec.grains.options.psf))
                    blobs(ii_d).psf = { permute(self.parameters.rec.grains.options.psf, [1 3 2]) };
                end

                if (strcmpi(self.parameters.acq(det_ind).type, 'topotomo'))
                    % This holds true even in case of single topotomo
                    % reconstruction: in that case the sampling size in z
                    % will be equal to 1
                    lambda_det(ii_d) = 1 / sampler.sampling_size(3);
                end

                geometries{ii_d} = algo_params.geometries;
                offsets{ii_d} = algo_params.offsets;
                proj_uv_size(ii_d, :) = algo_params.proj_uv_size;
            end

            algo = Gt6DBlobReconstructor(volume_size, blobs, proj_uv_size, ...
                'data_type', self.data_type, ...
                'geometries', geometries, ...
                'offsets', offsets, ...
                'volume_ss', self.volume_downscaling, ...
                'rspace_ss', self.rspace_super_sampling, ...
                'orientation_groups', or_groups, ...
                'lambda_det', lambda_det, ...
                varargin{:} );
        end

        function [algo, blobs_struct] = get6DAlgorithmMultiRegion(self, sampler, num_interp, blobs, varargin)
            % Forcing use_matrix_row_rescaling and
            % use_predicted_scatter_ints to remove problems due to overlaps
            self.use_matrix_row_rescaling = true;
            self.use_predicted_scatter_ints = max(self.use_predicted_scatter_ints, 1);

            blobs = reshape(blobs, [], 1);

            num_regions = numel(sampler);

            det_ind = self.det_index(1);

            % Build Empty volumes
            ref_gr = sampler(1).get_reference_grain();
            volume_size = self.get_volume_size(ref_gr, det_ind);
            vol_size_mm = volume_size .* self.parameters.recgeo(1).voxsize * self.volume_downscaling;

            blobs_selected = false(size(blobs));
            blobs_w_interp = zeros(size(blobs));
            num_interp_ors = zeros(num_regions, 1);

            tot_all_ors = 0;
            or_ranges = zeros(2, num_regions);

            fprintf('Computing blobs interpolation widths..')
            c = tic();
            for ii_o_reg = 1:num_regions
                % Deciding the type of num_interp type
                num_interp_ors(ii_o_reg) = self.compute_num_interp(num_interp, sampler(ii_o_reg), det_ind);

                blobs_w_interp_or = self.compute_blobs_w_interp(sampler(ii_o_reg), num_interp_ors(ii_o_reg), det_ind);

                local_sel = sampler(ii_o_reg).selected;

                gr = sampler(ii_o_reg).get_reference_grain();
                or_inc_pos = gr.proj(det_ind).global_pos(local_sel);
                blobs_selected(or_inc_pos) = true;
                % Here we remove the problem of different w-interp, by
                % simply selecting the highest
                blobs_w_interp(or_inc_pos) = max(blobs_w_interp(or_inc_pos), blobs_w_interp_or);

                num_ors_reg = sampler(ii_o_reg).get_total_num_orientations();
                or_ranges(:, ii_o_reg) = [1, num_ors_reg] + tot_all_ors;
                tot_all_ors = tot_all_ors + num_ors_reg;
            end
            fprintf('\b\b: Done in %g seconds.\n', toc(c));

            inc_to_sel = zeros(size(blobs_selected));
            inc_to_sel(blobs_selected) = 1:sum(blobs_selected);

            fprintf('Extracting blobs on detector..')
            c = tic();
            blobs = blobs(blobs_selected);
            fprintf('\b\b: Done in %g seconds.\n', toc(c));

            blobs_w_interp = blobs_w_interp(blobs_selected);

            ref_omegas = cat(1, blobs(:).bbwim);
            ref_omegas = sum(ref_omegas, 2) / 2 * gtAcqGetOmegaStep(self.parameters, det_ind);

            shape_funcs = cell(num_regions, 1);
            ors_allblobs = cell(num_regions, 1);
            ors_w_tab = cell(num_regions, 1);
            scatter_ints_avg = cell(num_regions, 1);
            scatter_ints_ors = cell(num_regions, 1);
            paddings_w = NaN(numel(blobs), 2);
            uv_tab = NaN(numel(blobs), 2, 2, tot_all_ors);

            fprintf('Computing projection limits for each orientation bbox:\n')
            for ii_o_reg = 1:num_regions
                gr = sampler(ii_o_reg).get_reference_grain();

                fprintf('%d) Grainid %d:\n', ii_o_reg, gr.id)

                local_sel = sampler(ii_o_reg).selected;

                or_inc_pos = gr.proj(det_ind).global_pos(local_sel);
                or_sel_pos = inc_to_sel(or_inc_pos);

                local_bls_w_interp = blobs_w_interp(or_sel_pos);
                local_bls = blobs(or_sel_pos);

                shape_funcs{ii_o_reg} = self.get_shape_functions(sampler(ii_o_reg), num_interp, local_bls_w_interp);

                [orients, orients_ss] = sampler(ii_o_reg).get_orientations();

                num_ospace_oversampling = prod(sampler(ii_o_reg).ss_factor);
                if (num_ospace_oversampling == 1)
                    ors_struct = [orients{:}];
                else
                    ors_struct = cat(4, orients_ss{:});
                    ors_struct = [ors_struct{:}];
                end
                % Array of structures that contain all the info relative to
                % each orientation, for each blob
                ors_allblobs{ii_o_reg} = [ors_struct(:).allblobs];

                fprintf(' - Dealing with paddings:\n')
                % Finding the projection limits in W, and taking care of the
                % 360 degrees wrapping
                switch (lower(self.shape_functions_type))
                    case {'none', 'n'}
                        fprintf('   * Computing W num-interp\n')
                        [w_tab_lims, ors_w_tab{ii_o_reg}, ref_ws] = self.get_w_tab_numinterp(gr, ors_allblobs{ii_o_reg}, ref_omegas(or_sel_pos), det_ind);
                    case {'w', 'nw'}
                        fprintf('   * Computing W shape-functions\n')
                        [w_tab_lims, ref_ws] = self.get_w_tab_shapefuncs(shape_funcs{ii_o_reg}, ref_omegas(or_sel_pos), det_ind);
                end

                fprintf('   * Computing W paddings\n')
                [blobs_w_lims_pad_or, projs_w_lims, paddings_w_or] = ...
                    self.get_blob_w_lims(w_tab_lims, local_bls, ref_ws, local_bls_w_interp);

                fprintf('   * Computing UV tab..')
                % <blobs x uv x [min, max] x orientations>
                c = tic();
                [uv_tab_or, ~] = self.get_uv_tab(gr, ors_allblobs{ii_o_reg}, vol_size_mm, det_ind);
                fprintf('\b\b: Done in %g seconds.\n', toc(c));

                sel_incl_indx = find(local_sel);

                local_bls_sizes = cat(1, local_bls(:).bbsize);
                local_bls_sizes = local_bls_sizes(:, 3) + sum(paddings_w_or, 2);
                local_bls_sizes = (local_bls_sizes - 1) ./ local_bls_w_interp + 1;

                [scatter_ints_avg{ii_o_reg}, scatter_ints_ors{ii_o_reg}] = self.get_scattering_intensities(gr, ors_allblobs{ii_o_reg}, det_ind);

                fprintf('------------------------+------------------------+--------------+---------+----------+------+----------+-------------\n')
                fprintf('Blob (n.)               | Blob w lims            | Projs w lims | Padding | n-interp | Size | Blob Int | Scatter Int\n')
                fprintf('------------------------+------------------------+--------------+---------+----------+------+----------+-------------\n')
                for ii_b = 1:numel(local_bls)
                    fprintf('%03d (%03d -> %03d -> %03d) | %4d:%4d -> %4d:%4d |    %4d:%4d |   %2d:%2d |       %2d | 1:%2d | %8s | %s [%s %s]\n', ...
                        ii_b, sel_incl_indx(ii_b), or_inc_pos(ii_b), ...
                        or_sel_pos(ii_b), local_bls(ii_b).bbwim, ...
                        blobs_w_lims_pad_or(ii_b, :),...
                        round(projs_w_lims(ii_b, :)), ...
                        paddings_w_or(ii_b, :), local_bls_w_interp(ii_b), ...
                        local_bls_sizes(ii_b), ...
                        sprintf('%3.3g', local_bls(ii_b).intensity), ...
                        sprintf('%3.3g', scatter_ints_avg{ii_o_reg}(ii_b)), ...
                        sprintf('%3.3g', min(scatter_ints_ors{ii_o_reg}(ii_b, :))), ...
                        sprintf('%3.3g', max(scatter_ints_ors{ii_o_reg}(ii_b, :))) );
                end
                fprintf('------------------------+------------------------+--------------+---------+----------+------+----------+-------------\n')

                % Updating paddings required by this orientation box (since
                % some indices might repeat, we do it in a for loop)
                for ii_b = 1:numel(or_sel_pos)
                    paddings_w(or_sel_pos(ii_b), :) = max(paddings_w(or_sel_pos(ii_b), :), paddings_w_or(ii_b, :));

                    uv_tab(or_sel_pos(ii_b), :, :, or_ranges(1, ii_o_reg):or_ranges(2, ii_o_reg)) = uv_tab_or(ii_b, :, :, :);
                end
                fprintf('\n')
            end

            lims_blobs_w_orig = cat(1, blobs(:).bbwim);
            blobs_w_lims_pad = [ ...
                lims_blobs_w_orig(:, 1) - paddings_w(:, 1), ...
                lims_blobs_w_orig(:, 2) + paddings_w(:, 2), ...
                ];

            fprintf('Computing UV paddings..')
            [blobs_uv_lims, blob_padding_uv, proj_uv_size, projs_uv_lims, shifts_uv] = self.get_blob_uv_lims(blobs, uv_tab);
            fprintf('\b\b => blob sizes UV: [%d, %d] -> [%d, %d], projection sizes UV: [%d, %d]\n\n', ...
                blobs(1).bbsize(1:2), blobs(1).bbsize(1:2) + sum(blob_padding_uv, 1), proj_uv_size)

            bls = self.pad_blobs(blobs, paddings_w, blobs_w_interp, blob_padding_uv);

            blobs_sizes_w = zeros(1, numel(bls));
            % They were renormalized to 1 in pad_blobs
            for ii_b = 1:numel(bls)
                bls{ii_b} = bls{ii_b} * blobs(ii_b).intensity;

                blobs_sizes_w(ii_b) = size(bls{ii_b}, 2);
            end

            proj_props = struct( ...
                'geom', cell(num_regions, 1), ...
                'offs', cell(num_regions, 1), ...
                'geom_ss', cell(num_regions, 1), ...
                'offs_ss', cell(num_regions, 1) );

            fprintf('Completing projection geometry determination:\n');
            for ii_o_reg = 1:num_regions
                gr = sampler(ii_o_reg).get_reference_grain();

                local_sel = sampler(ii_o_reg).selected;

                if (~isfield(gr.proj(det_ind), 'global_pos') ...
                        || isempty(gr.proj(det_ind).global_pos))
                    gr.proj(det_ind).global_pos = zeros(size(sampler(ii_o_reg).ondet));
                    gr.proj(det_ind).global_pos = 1:numel(local_sel);
                end
                or_inc_pos = gr.proj(det_ind).global_pos(local_sel);
                or_sel_pos = inc_to_sel(or_inc_pos);

                local_bls_w_interp = blobs_w_interp(or_sel_pos);
                local_bls_w_lims_pad = blobs_w_lims_pad(or_sel_pos, :);
                % Dimensions are: <blobs x uv x [min, max] x orientations>
                local_projs_uv_lims = projs_uv_lims(or_sel_pos, :, :, or_ranges(1, ii_o_reg):or_ranges(2, ii_o_reg));
                % Dimensions are: <blobs x uv x orientations>
                local_shifts_uv = shifts_uv(or_sel_pos, :, or_ranges(1, ii_o_reg):or_ranges(2, ii_o_reg));

                fprintf(' - Computing ASTRA proj geometry and blobs <-> sinograms coefficients..')
                c = tic();
                % ASTRA Geometry
                switch (lower(self.shape_functions_type))
                    case 'none'
                        geometries = self.compute_proj_geom_numinterp( ...
                            local_projs_uv_lims, gr, ors_allblobs{ii_o_reg}, det_ind);

                        offsets = self.compute_proj_coeffs_numinterp(...
                            local_bls_w_lims_pad, ors_w_tab{ii_o_reg}, local_bls_w_interp, local_shifts_uv);
                    case 'n'
                        error('Gt6DReconstructionAlgorithmFactory:wrong_argument', ...
                            'Projection geometry with only N-shape-functions, is not implemented, yet!')
                    case 'w'
                        [geometries, offsets] = self.compute_proj_geom_shapefuncs_w( ...
                            local_projs_uv_lims, gr, ors_allblobs{ii_o_reg}, ...
                            local_bls_w_lims_pad, shape_funcs{ii_o_reg}, local_bls_w_interp, local_shifts_uv, det_ind);
                    case 'nw'
                        [geometries, offsets] = self.compute_proj_geom_shapefuncs_nw( ...
                            local_projs_uv_lims, gr, ors_allblobs{ii_o_reg}, ...
                            local_bls_w_lims_pad, shape_funcs{ii_o_reg}, local_bls_w_interp, local_shifts_uv, det_ind);
                end
                fprintf('\b\b: Done in %g seconds.\n', toc(c));

                fprintf(' - Rescaling projection coefficients, based on scattering intensity..')
                c = tic();
                for ii_o = 1:numel(offsets)
                    for ii_b = 1:numel(scatter_ints_avg{ii_o_reg})
                        inds = offsets{ii_o}.blob_offsets == ii_b;
                        offsets{ii_o}.proj_coeffs(inds) = offsets{ii_o}.proj_coeffs(inds) .* scatter_ints_ors{ii_o_reg}(ii_b, ii_o);
                    end
                end
                fprintf('\b\b: Done in %g seconds.\n', toc(c));

                fprintf(' - Redirecting blobs positions..')
                c = tic();
                for ii_o = 1:numel(offsets)
                    offsets{ii_o}.blob_offsets = or_sel_pos(offsets{ii_o}.blob_offsets);
                end
                fprintf('\b\b: Done in %g seconds.\n', toc(c));

                % Taking care of the orientation-space super-sampling
                tot_orient = sampler(ii_o_reg).get_total_num_orientations();
                if (num_ospace_oversampling > 1)
                    geometries = reshape(geometries, num_ospace_oversampling, []);
                    offsets = reshape(offsets, num_ospace_oversampling, []);

                    geometries_ss = cell(tot_orient, 1);
                    offsets_ss = cell(tot_orient, 1);
                    for ii_or = 1:tot_orient
                        geometries_ss{ii_or} = cat(1, geometries{:, ii_or});

                        % merging all the orientation-space super-sampling
                        % projection coefficients into one structure per
                        % orientation. This is needed in order to fwd/bwd-project
                        % all the different sub-orientations together.
                        offsets_ss{ii_or} = self.merge_offsets(offsets(:, ii_or));
                    end

                    geometries = geometries_ss;
                    offsets = offsets_ss;
                end

                proj_props(ii_o_reg).geom = geometries;
                proj_props(ii_o_reg).offs = offsets;
            end

            geometries = cat(1, proj_props(:).geom);
            offsets = cat(1, proj_props(:).offs);

            % Handling PSFs (very simplistic ATM)
            if (isfield(self.parameters.rec.grains.options, 'psf') ...
                    && ~isempty(self.parameters.rec.grains.options.psf))
                fprintf('\nGlobal PSF used for all orientation-space regions\n\n')
                psf = { permute(self.parameters.rec.grains.options.psf, [1 3 2]) };
            else
                fprintf('\nNo PSF used for any orientation-space region\n\n')
                psf = {};
            end

            p = self.parameters;
            pixel_size = [p.detgeo(det_ind).pixelsizeu, p.detgeo(det_ind).pixelsizev];
            pixel_size_ratio = pixel_size ./ mean(p.recgeo(self.det_index(1)).voxsize);
            proj_det_ss = mean(pixel_size_ratio);

            blobs_struct = gt6DRecBlobsDefinition();
            blobs_struct(1).data = bls;
            blobs_struct(1).bbuims = squeeze(blobs_uv_lims(:, 1, :));
            blobs_struct(1).bbvims = squeeze(blobs_uv_lims(:, 2, :));
            blobs_struct(1).bbwims = blobs_w_lims_pad;
            blobs_struct(1).size_uv = [size(bls{ii_b}, 1), size(bls{ii_b}, 3)];
            blobs_struct(1).sizes_w = blobs_sizes_w;
            blobs_struct(1).det_index = det_ind;
            blobs_struct(1).det_ss = proj_det_ss;
            blobs_struct(1).psf = psf;

            or_groups = self.get_orientation_groups(sampler);

            algo = Gt6DBlobReconstructor(volume_size, blobs_struct, proj_uv_size,...
                'data_type', self.data_type, ...
                'geometries', {geometries}, ...
                'offsets', {offsets}, ...
                'volume_ss', self.volume_downscaling, ...
                'orientation_groups', or_groups, ...
                varargin{:} );
        end
    end

    methods (Access = protected)
        function shape_funcs = get_shape_functions(self, sampler, num_interp, blobs_w_interp, det_ind)
            if (~exist('det_ind', 'var') || isempty(det_ind))
                if (isempty(self.det_index))
                    det_ind = 1;
                else
                    det_ind = self.det_index(1);
                end
            end

            switch (lower(self.shape_functions_type))
                case 'none'
                    shape_funcs = {};
                case 'n'
%                     shape_funcs = gtDefShapeFunctionsCreateNW(sampler);
                    shape_funcs = gtDefShapeFunctionsFwdProj(sampler, ...
                        'shape_function_type', 'nw', ...
                        'num_interp', num_interp, ...
                        'blobs_w_interp', blobs_w_interp, ...
                        'det_ind', det_ind);
                    shape_funcs = gtDefShapeFunctionsNW2N(shape_funcs);
                case 'w'
%                     shape_funcs = gtDefShapeFunctionsCreateW(sampler);
                    shape_funcs = gtDefShapeFunctionsFwdProj(sampler,  ...
                        'shape_function_type', 'w', ...
                        'num_interp', num_interp, ...
                        'blobs_w_interp', blobs_w_interp, ...
                        'det_ind', det_ind);
                case 'nw'
%                     shape_funcs = gtDefShapeFunctionsCreateNW(sampler);
                    shape_funcs = gtDefShapeFunctionsFwdProj(sampler, ...
                        'shape_function_type', 'nw', ...
                        'num_interp', num_interp, ...
                        'blobs_w_interp', blobs_w_interp, ...
                        'det_ind', det_ind);
                otherwise
                    error('Gt6DReconstructionAlgorithmFactory:wrong_argument', ...
                        'Shape functions of type: %s, not supported yet!', ...
                        self.shape_functions_type)
            end
        end

        function [shared, parent_corresp_shared_blobs, shared_pos] = get_shared_bls(self, selected_pos_in_included, sampler, ii_o_reg)
            curr_or = sampler(ii_o_reg).get_reference_grain();
            curr_sel = sampler(ii_o_reg).selected;
            % This contains the position of the shared blobs in
            % parent's numbering
            shared_pos = curr_or.proj(self.det_index).shared_bls_with_parent(curr_sel);
            % This indicates the ones that are actually shared among
            % the twin's blobs
            shared = shared_pos > 0;

            parent_corresp_shared_blobs = selected_pos_in_included(shared_pos(shared));
        end

        function algo_params = get_algorithm_params(self, sampler, blobs_w_interp, shape_funcs, vol_size, det_ind)
            is_topotomo = strcmpi(self.parameters.acq(det_ind).type, 'topotomo');

            ref_gr = sampler.get_reference_grain();

            [orients, orients_ss] = sampler.get_orientations(det_ind);

            num_ospace_oversampling = prod(sampler.ss_factor);
            if (num_ospace_oversampling == 1)
                ors_struct = [orients{:}];
            else
                ors_struct = cat(4, orients_ss{:});
                ors_struct = [ors_struct{:}];
            end
            % Array of structures that contain all the info relative to
            % each orientation, for each blob
            ors_allblobs = [ors_struct(:).allblobs];

            ond = ref_gr.proj(det_ind).ondet;
            inc = ref_gr.proj(det_ind).included;
            sel = ref_gr.proj(det_ind).selected;

            fprintf(' - Extracting blobs: ')
            c = tic();
            bl = ref_gr.proj(det_ind).bl(sel);
            fprintf('\b\b (%f s).\n', toc(c));
            num_blobs = numel(bl);
            sel_incl_indx = find(sel);

            fprintf(' - Dealing with blobs/scattering intensities:\n')
            if (self.use_predicted_scatter_ints)
                fprintf('   * Computing families'' scattering intensities\n')
                [scatter_ints_avg, scatter_ints_ors] = self.get_scattering_intensities(ref_gr, ors_allblobs, det_ind);
            else
                fprintf('   * Using blobs'' intensities\n')
                scatter_ints_avg = [bl(:).intensity]';
                scatter_ints_avg = scatter_ints_avg ./ mean([bl(:).intensity]);
                scatter_ints_ors = scatter_ints_avg(:, ones(1, numel(ors_allblobs)));
            end

            fprintf(' - Dealing with paddings:\n')
            if (~is_topotomo)
                ref_omegas = ref_gr.allblobs(det_ind).omega(ond(inc(sel)));

                % Finding the projection limits in W, and taking care of the
                % 360 degrees wrapping
                switch (lower(self.shape_functions_type))
                    case {'none', 'n'}
                        fprintf('   * Computing W num-interp\n')
                        [w_tab_lims, w_tab, ref_ws] = self.get_w_tab_numinterp(ref_gr, ors_allblobs, ref_omegas, det_ind);
                    case {'w', 'nw'}
                        fprintf('   * Computing W shape-functions\n')
                        [w_tab_lims, ref_ws] = self.get_w_tab_shapefuncs(shape_funcs, ref_omegas, det_ind);
                end

                fprintf('   * Computing W paddings\n')
                [blobs_w_lims_pad, projs_w_lims, paddings] = ...
                    self.get_blob_w_lims(w_tab_lims, bl, ref_ws, blobs_w_interp);
            else
                fprintf('   * Computing P num-interp\n')
                [w_tab_lims, w_tab] = self.get_p_tab_numinterp(ref_gr, ors_allblobs, det_ind);

                fprintf('   * Computing P paddings\n')
                [blobs_w_lims_pad, projs_w_lims, paddings] = ...
                    self.get_blob_p_lims(w_tab_lims, bl, blobs_w_interp);
            end

            fprintf('   * Computing UV paddings..')
            [uv_tab, ~] = self.get_uv_tab(ref_gr, ors_allblobs, vol_size, det_ind);
            [blobs_uv_lims, blob_padding_uv, proj_uv_size, projs_uv_lims, shifts_uv] = self.get_blob_uv_lims(bl, uv_tab);

            fprintf('\b\b => blob sizes UV: [%d, %d] -> [%d, %d], projection sizes UV: [%d, %d]\n', ...
                bl(1).bbsize(1:2), bl(1).bbsize(1:2) + sum(blob_padding_uv, 1), proj_uv_size)

            blobs = self.pad_blobs(bl, paddings, blobs_w_interp, blob_padding_uv);
            if (self.use_matrix_row_rescaling)
                renorm_int = [bl(:).intensity]';
            elseif (mean([bl(:).intensity]) < 10)
                % If dealing with old style blobs
                pix_size = self.parameters.recgeo(self.det_index(1)).voxsize * self.volume_downscaling;
                renorm_int = prod(vol_size ./ pix_size) / 10;
                renorm_int = ones(num_blobs, 1) * renorm_int;
            else
                renorm_int = [bl(:).intensity]' ./ scatter_ints_avg;
            end

            blobs_sizes_w = zeros(1, num_blobs);
            blobs_ints = zeros(1, num_blobs);
            % They were renormalized to 1 in pad_blobs
            for ii_b = 1:numel(blobs)
                blobs{ii_b} = blobs{ii_b} .* renorm_int(ii_b);
                blobs_ints(ii_b) = sum(blobs{ii_b}(:));

                blobs_sizes_w(ii_b) = size(blobs{ii_b}, 2);
            end
            if (is_topotomo)
                blobs_orig_w_lims = cat(1, bl(:).bbpim);
            else
                blobs_orig_w_lims = cat(1, bl(:).bbwim);
            end

            % Let's restrict blobs to the reached regions
            fprintf('----------+------------------------+--------------+---------+----------+------+----------+-------------\n')
            fprintf('Blob (n.) | Blob w lims            | Projs w lims | Padding | n-interp | Size | Blob Int | Scatter Int\n')
            fprintf('----------+------------------------+--------------+---------+----------+------+----------+-------------\n')
            for ii_b = 1:num_blobs
                try
                fprintf('%03d (%03d) | %4d:%4d -> %4d:%4d |    %4d:%4d |   %2d:%2d |       %2d | 1:%2d | %8s | %s [%s %s]\n', ...
                    ii_b, sel_incl_indx(ii_b), ...
                    blobs_orig_w_lims(ii_b, :), blobs_w_lims_pad(ii_b, :),...
                    round(projs_w_lims(ii_b, :)), ...
                    paddings(ii_b, :), blobs_w_interp(ii_b), ...
                    blobs_sizes_w(ii_b), ...
                    sprintf('%3.3g', gtMathsSumNDVol(blobs{ii_b})), ...
                    sprintf('%3.3g', blobs_ints(ii_b)), ...
                    sprintf('%3.3g', min(scatter_ints_ors(ii_b, :))), ...
                    sprintf('%3.3g', max(scatter_ints_ors(ii_b, :))) );
                catch mexc
                    gtPrintException(mexc)
                end
            end
            fprintf('----------+------------------------+--------------+---------+----------+------+----------+-------------\n')

            fprintf(' - Computing ASTRA proj geometry and blobs <-> sinograms coefficients: ')
            c = tic();
            % ASTRA Geometry
            switch (lower(self.shape_functions_type))
                case 'none'
                    geometries = self.compute_proj_geom_numinterp( ...
                        projs_uv_lims, ref_gr, ors_allblobs, det_ind);

                    offsets = self.compute_proj_coeffs_numinterp(...
                        blobs_w_lims_pad, w_tab, blobs_w_interp, shifts_uv);
                case 'n'
                    error('Gt6DReconstructionAlgorithmFactory:wrong_argument', ...
                        'Projection geometry with only N-shape-functions, is not implemented, yet!')
                case 'w'
                    [geometries, offsets] = self.compute_proj_geom_shapefuncs_w( ...
                        projs_uv_lims, ref_gr, ors_allblobs, blobs_w_lims_pad, ...
                        shape_funcs, blobs_w_interp, shifts_uv, det_ind);
                case 'nw'
                    [geometries, offsets] = self.compute_proj_geom_shapefuncs_nw( ...
                        projs_uv_lims, ref_gr, ors_allblobs, blobs_w_lims_pad, ...
                        shape_funcs, blobs_w_interp, shifts_uv, det_ind);
            end
            fprintf('Done (%f s).\n', toc(c));

            fprintf(' - Rescaling projection coefficients, based on scattering intensity: ')
            c = tic();
            if (self.use_matrix_row_rescaling)
                for ii_o = 1:numel(offsets)
                    for ii_b = 1:numel(scatter_ints_avg)
                        inds = offsets{ii_o}.blob_offsets == ii_b;
                        offsets{ii_o}.proj_coeffs(inds) = offsets{ii_o}.proj_coeffs(inds) .* scatter_ints_ors(ii_b, ii_o);
                    end
                end
            end
            fprintf('Done (%f s).\n', toc(c));

            % Taking care of the orientation-space super-sampling
            tot_orient = sampler.get_total_num_orientations();
            if (num_ospace_oversampling > 1)
                geometries = reshape(geometries, num_ospace_oversampling, []);
                offsets = reshape(offsets, num_ospace_oversampling, []);

                geometries_ss = cell(tot_orient, 1);
                offsets_ss = cell(tot_orient, 1);
                for ii_or = 1:tot_orient
                    geometries_ss{ii_or} = cat(1, geometries{:, ii_or});

                    % merging all the orientation-space super-sampling
                    % projection coefficients into one structure per
                    % orientation. This is needed in order to fwd/bwd-project
                    % all the different sub-orientations together.
                    offsets_ss{ii_or} = self.merge_offsets(offsets(:, ii_or));
                end

                geometries = geometries_ss;
                offsets = offsets_ss;
            end

            p = self.parameters;
            pixel_size = [p.detgeo(det_ind).pixelsizeu, p.detgeo(det_ind).pixelsizev];
            pixel_size_ratio = pixel_size ./ mean(p.recgeo(self.det_index(1)).voxsize);
            proj_det_ss = mean(pixel_size_ratio);

            if (isfield(bl, 'psf') && ~isempty(bl.psf))
                psf = arrayfun(@(x){ permute(x.psf, [1 3 2]) }, bl);
            elseif (isfield(ref_gr.proj(det_ind), 'psf') && ~isempty(ref_gr.proj(det_ind).psf))
                psf = { permute(ref_gr.proj(det_ind).psf, [1 3 2]) };
            else
                psf = {};
            end

            blobs_struct = gt6DRecBlobsDefinition();
            blobs_struct.data = blobs;
            blobs_struct.bbuims = squeeze(blobs_uv_lims(:, 1, :));
            blobs_struct.bbvims = squeeze(blobs_uv_lims(:, 2, :));
            blobs_struct.bbwims = blobs_w_lims_pad;
            blobs_struct.size_uv = [size(blobs{ii_b}, 1), size(blobs{ii_b}, 3)];
            blobs_struct.sizes_w = blobs_sizes_w;
            blobs_struct.det_index = det_ind;
            blobs_struct.det_ss = proj_det_ss;
            blobs_struct.psf = psf;

            algo_params = struct( ...
                'blobs', {blobs_struct}, ...
                'geometries', {geometries}, ...
                'offsets', {offsets}, ...
                'proj_uv_size', proj_uv_size, ...
                'blobs_w_lims', {blobs_w_lims_pad} );
        end

        function volume_size = get_volume_size(self, ref_gr, det_ind)
            proj = ref_gr.proj(det_ind);

            spacing = mean([proj.vol_size_y, proj.vol_size_x, proj.vol_size_z]) * (self.rspace_oversize - 1);
            volume_size = ceil([proj.vol_size_y, proj.vol_size_x, proj.vol_size_z] + spacing);

            if (self.volume_downscaling > 1)
                volume_size = ceil(volume_size / self.volume_downscaling);
            end
        end

        function num_interp = compute_num_interp(~, num_interp, sampler, det_ind)
            if (num_interp < 0)
                num_interp = abs(num_interp) * sampler.minimum_num_interp(det_ind);
            elseif (~num_interp)
                num_interp = sampler.suggest_num_interp(det_ind);
            end
        end

        function blobs_w_interp = compute_blobs_w_interp(self, sampler, num_interp, det_ind)
            is_topotomo = strcmpi(self.parameters.acq(det_ind).type, 'topotomo');

            ref_or = sampler.get_reference_grain();
            ond = ref_or.proj(det_ind).ondet;
            inc = ref_or.proj(det_ind).included;
            sel = ref_or.proj(det_ind).selected;
            ab_sel = ond(inc(sel));

            if (strcmpi(self.num_interp_type, 'lorentz') && ~is_topotomo)
                etas = ref_or.allblobs(det_ind).eta(ab_sel);
                lorentz_factor = 1 ./ abs(sind(etas));
                % Let's try to not penalize rounding errors while being a
                % bit tollerant to small deviations from integer numbers
                % (rounding in [-0.75 0.25), instead of [-0.5 0.5) ), but
                % still trying to avoid gaps in the projections
                blobs_w_interp = ceil(num_interp .* (lorentz_factor - 0.25));
            else
                blobs_w_interp = ceil(num_interp(ones(numel(ab_sel), 1), 1));
            end
        end

        function pad_interp_blobs = pad_blobs(self, blobs, paddings_w, slices_interp, padding_uv)
            num_blobs = numel(blobs);
            pad_interp_blobs = cell(num_blobs, 1);

            blob_original_sizes = cat(1, blobs(:).bbsize);

            % nums_pad_slices : number of Ws in the padded blob
            % nums_interp_slices : number of groupped slices in the exit
            %                      blobs
            nums_pad_slices = blob_original_sizes(:, 3) + sum(paddings_w, 2);
            nums_interp_slices = (nums_pad_slices -1) ./ slices_interp + 1;

            blobs_uv_sizes = bsxfun(@plus, blob_original_sizes(:, 1:2), sum(padding_uv, 1));
            blobs_sizes = [ ...
                blobs_uv_sizes(:, 1), nums_interp_slices, blobs_uv_sizes(:, 2)];

            fprintf(' - Creating padded and interpolated blobs: ')
            c = tic();
            for ii = 1:num_blobs
                num_chars = fprintf('%03d/%03d', ii, num_blobs);

                blob = blobs(ii).intm;
                blob(blob < 0) = 0;
                nan_inds = isnan(blob);
                if (any(nan_inds(:)))
                    warning('Gt6DReconstructionAlgorithmFactory:bad_blob', ...
                        'Blob %d contains NaN! Setting them to 0', ii)
                end
                blob(nan_inds) = 0;
                blob = blob ./ gtMathsSumNDVol(blob);

                pad_interp_blobs{ii} = zeros(blobs_sizes(ii, :), self.data_type);

                slice_inds_u = (1:blob_original_sizes(ii, 1)) + padding_uv(1, 1);
                slice_inds_v = (1:blob_original_sizes(ii, 2)) + padding_uv(1, 2);

                % Let's build the interpolated slices
                for w_pad = 1:nums_interp_slices(ii)
                    coeffs = abs( (slices_interp(ii) * (w_pad-1) +1) - (1:nums_pad_slices(ii)));
                    coeffs = 1 - coeffs / slices_interp(ii);

                    % these are the padded (non interpolated) blob slices
                    % that relate to the current ii interpolated slice
                    blob_inds = find(coeffs > 0);
                    % The corresponding coefficients
                    coeffs = coeffs(blob_inds);

                    % The indices reported in the original blob (without
                    % padding)
                    blob_inds = blob_inds - paddings_w(ii, 1);
                    safe_blob_inds_pos = (blob_inds <= blob_original_sizes(ii, 3)) & (blob_inds >= 1);
                    safe_blob_inds = blob_inds(safe_blob_inds_pos);
                    safe_coeffs = coeffs(safe_blob_inds_pos);

                    % We now produce the slice
                    slice = zeros(blob_original_sizes(ii, 1:2), self.data_type);
                    for w = 1:numel(safe_coeffs)
                        slice = slice + safe_coeffs(w) .* blob(:, :, safe_blob_inds(w));
                    end
                    pad_interp_blobs{ii}(slice_inds_u, w_pad, slice_inds_v) = slice;
                end

                fprintf(repmat('\b', [1 num_chars]));
            end
            fprintf('\b\b (%f s).\n', toc(c));
        end

        function [w_tab_lims, w_tab, ref_w] = get_w_tab_numinterp(self, ref_gr, or_allblobs, ref_ws, det_ind)
            ond = ref_gr.proj(det_ind).ondet;
            inc = ref_gr.proj(det_ind).included;
            sel = ref_gr.proj(det_ind).selected;

            ab_sel = ond(inc(sel));

            omega_step = gtAcqGetOmegaStep(self.parameters, det_ind);
            % Tabular of omegas: horizontal direction the n geometries,
            % vertical direction the omegas relative to each blob
            w_tab = [or_allblobs(:).omega];
            % Images numbers
            w_tab = w_tab(ab_sel, :) / omega_step;

            ref_w = ref_ws / omega_step;

            w_tab = gtGrainAnglesTabularFix360deg(w_tab, ref_w, self.parameters, det_ind);
            w_tab_lims = cat(3, floor(w_tab), ceil(w_tab));
        end

        function [p_tab_lims, p_tab] = get_p_tab_numinterp(self, ref_gr, or_allblobs, det_ind)
            ond = ref_gr.proj(det_ind).ondet;
            inc = ref_gr.proj(det_ind).included;
            sel = ref_gr.proj(det_ind).selected;

            ab_sel = ond(inc(sel));

            ph_step = gtAcqGetPhiStep(self.parameters, det_ind);
            % Tabular of phis: horizontal direction the n geometries,
            % vertical direction the omegas relative to each blob
            p_tab = [or_allblobs(:).phi];
            % Images numbers
            p_tab = p_tab(ab_sel, :) / ph_step;

            p_tab_lims = cat(3, floor(p_tab), ceil(p_tab));
        end

        function [w_tab_lims, ref_w] = get_w_tab_shapefuncs(self, shape_funcs, ref_ws, det_ind)
            omega_step = gtAcqGetOmegaStep(self.parameters, det_ind);

            num_ors = numel(shape_funcs);
            num_blobs = numel(shape_funcs{1});
            w_tab_lims = zeros(num_blobs, num_ors, 2);
            % Tabular of omegas: horizontal direction the n geometries,
            % vertical direction the omegas relative to each blob
            for ii_or = 1:num_ors
                w_tab_lims(:, ii_or, :) = cat(1, shape_funcs{ii_or}(:).bbwim);
            end

            ref_w = ref_ws / omega_step;

            w_tab_lims(:, :, 1) = gtGrainAnglesTabularFix360deg(w_tab_lims(:, :, 1), ref_w, self.parameters, det_ind);
            w_tab_lims(:, :, 2) = gtGrainAnglesTabularFix360deg(w_tab_lims(:, :, 2), ref_w, self.parameters, det_ind);
        end

        function [uv_tab, uv_tab_lims] = get_uv_tab(self, ref_gr, ors_allblobs, vol_size, det_ind)
            ond = ref_gr.proj(det_ind).ondet;
            inc = ref_gr.proj(det_ind).included;
            sel = ref_gr.proj(det_ind).selected;

            ab_sel = ond(inc(sel));

            num_ors = numel(ors_allblobs);
            num_blobs = numel(ab_sel);

            detgeo = self.parameters.detgeo(det_ind);
            detref_uv = [detgeo.detrefu, detgeo.detrefv]';
            labgeo = self.parameters.labgeo;
            samgeo = self.parameters.samgeo;

            do_correct_sam_shifts = isfield(self.parameters.acq, 'correct_sample_shifts') ...
                && (self.parameters.acq.correct_sample_shifts);
            if (do_correct_sam_shifts)
                om = zeros(numel(ab_sel), num_ors);
                for ii_o = 1:num_ors
                    om(:, ii_o) = ors_allblobs(ii_o).omega(ab_sel);
                end
                [~, shifts_sam] = gtMatchGetSampleShifts(self.parameters, om(:));
            end

            % Dimensions are: <blobs x uv x edges x orientations>
            uv_tab = zeros(num_blobs, 2, 8, num_ors);
            for ii_o = 1:num_ors
                dveclab_t = ors_allblobs(ii_o).dvec(ab_sel, :)';

                rot_s2l = ors_allblobs(ii_o).rot_s2l(:, :, ab_sel);
                gcsam_v = ref_gr.center(ones(1, size(rot_s2l, 3)), :);

                if (do_correct_sam_shifts)
                    % Adding sample shifts
                    shifts_sam_ii_o = shifts_sam( ((ii_o - 1) * numel(ab_sel) +1) : (ii_o * numel(ab_sel)), :);
                    gcsam_v = gcsam_v + shifts_sam_ii_o;
                end

                ii_e = 1;

                for dx = -0.5:0.5
                    for dy = -0.5:0.5
                        for dz = -0.5:0.5
                            edge_sam_v = bsxfun(@plus, gcsam_v, [dx, dy, dz] .* vol_size);
                            edge_lab_v_t = gtGeoSam2Lab(edge_sam_v, rot_s2l, labgeo, samgeo, false)';

                            uv_tab(:, :, ii_e, ii_o) = gtFedPredictUVMultiple([], dveclab_t, edge_lab_v_t, ...
                                detgeo.detrefpos', detgeo.detnorm', detgeo.Qdet, ...
                                detref_uv)';

                            ii_e = ii_e + 1;
                        end
                    end
                end
            end
            % Dimensions are: <blobs x uv x [min, max] x orientations>
            uv_tab = cat(3, floor(min(uv_tab, [], 3)), ceil(max(uv_tab, [], 3)));

            % Dimensions are: <blobs x uv x [min, max]>
            uv_tab_lims = cat(3, ...
                min(uv_tab(:, :, 1, :), [], 4), ...
                max(uv_tab(:, :, 2, :), [], 4));
        end

        function [blobs_uv_lims, padding_uv, proj_uv_size_max, projs_uv_lims, shifts_uv] = get_blob_uv_lims(~, bl, uv_tab)
            bls_bbuim = cat(1, bl.bbuim);
            bls_bbvim = cat(1, bl.bbvim);

            % Dimensions are: <blobs x uv x [min, max]>
            blobs_uv_lims_orig = cat(2, reshape(bls_bbuim, [], 1, 2), reshape(bls_bbvim, [], 1, 2));

            % Here we compute the maximum projection size among all
            % orientations on each blob, and then use it to compute
            % paddings for the projections, to have equal projections size
            % for each reflections of each orientation.
            % We are assuming that the center of the [min, max] range is
            % approximately the center of the projection.
            % Otherwise, some more complicated heuristic involving the
            % uv_tab_limits should be involved

            % Dimensions are: <blobs x uv x orientations>
            proj_uv_sizes = squeeze(uv_tab(:, :, 2, :) - uv_tab(:, :, 1, :) + 1);
            % Dimensions are: <blobs x uv x 1>
            proj_uv_size_bl = max(proj_uv_sizes, [], 3);
            % Dimensions are: <1 x uv x 1>
            proj_uv_size_max = max(proj_uv_size_bl, [], 1);

            proj_uv_size_max = proj_uv_size_max + mod(4 - mod(proj_uv_size_max, 4), 4);

            diffs_proj_uv_sizes = bsxfun(@minus, proj_uv_size_max, proj_uv_sizes) / 2;
            paddings_uv_sizes_lower = ceil(diffs_proj_uv_sizes);
            paddings_uv_sizes_upper = floor(diffs_proj_uv_sizes);

            num_ors = size(uv_tab, 4);
            num_blobs = size(uv_tab, 1);

            % Dimensions are: <blobs x uv x [min, max] x orientations>
            projs_uv_lims = cat(3, ...
                uv_tab(:, :, 1, :) - reshape(paddings_uv_sizes_lower, num_blobs, 2, 1, num_ors), ...
                uv_tab(:, :, 2, :) + reshape(paddings_uv_sizes_upper, num_blobs, 2, 1, num_ors) );

            % Dimensions are: <blobs x uv x [min, max]>
            all_projs_uv_lims = cat(3, ...
                min(projs_uv_lims(:, :, 1, :), [], 4), ...
                max(projs_uv_lims(:, :, 2, :), [], 4) );

            % We now compute the required blob paddings, on the basis of
            % the new computed projection sizes
            paddings_bl_uv_lower = max(blobs_uv_lims_orig(:, :, 1) - all_projs_uv_lims(:, :, 1), 0);
            paddings_bl_uv_upper = max(all_projs_uv_lims(:, :, 2) - blobs_uv_lims_orig(:, :, 2), 0);

            % Dimensions are: <[min, max] x uv>
            padding_uv = cat(1, max(paddings_bl_uv_lower, [], 1), max(paddings_bl_uv_upper, [], 1));

            % Dimensions are: <blobs x uv x [min, max]>
            blobs_uv_lims = cat(3, ...
                bsxfun(@minus, blobs_uv_lims_orig(:, :, 1), padding_uv(1, :)), ...
                bsxfun(@plus, blobs_uv_lims_orig(:, :, 2), padding_uv(2, :)) );

            % Dimensions are: <blobs x uv x orientations>
            shifts_uv = bsxfun(@minus, squeeze(projs_uv_lims(:, :, 1, :)), blobs_uv_lims(:, :, 1));
        end

        function [lims_blobs_w_pad, lims_projs_w, paddings] = get_blob_w_lims(self, w_tab_lims, bl, ref_int_ws, slices_interp)
        % Finding limits for the shape functions

            % Let's add extreme values
            lims_projs_w = [ ...
                min([w_tab_lims(:, :, 1), w_tab_lims(:, :, 2)], [], 2), ...
                max([w_tab_lims(:, :, 1), w_tab_lims(:, :, 2)], [], 2) ];

            lims_blobs_w_orig = vertcat( bl(:).bbwim );

            % Let's treat those blobs at the w edge 360->0
            % (from the blobs perspective)
            lims_blobs_w_orig = gtGrainAnglesTabularFix360deg(lims_blobs_w_orig, ref_int_ws, self.parameters);

            paddings(:, 1) = max(lims_blobs_w_orig(:, 1) - lims_projs_w(:, 1), 0);
            paddings(:, 2) = max(lims_projs_w(:, 2) - lims_blobs_w_orig(:, 2), 0);

            % Adding the padding
            lims_blobs_w_pad = [ ...
                lims_blobs_w_orig(:, 1) - paddings(:, 1), ...
                lims_blobs_w_orig(:, 2) + paddings(:, 2) ];

            % computing the needed extension for the interpolation
            lims_blobs_w_pad = [ ...
                floor(lims_blobs_w_pad(:, 1) ./ slices_interp) .* slices_interp, ...
                ceil(lims_blobs_w_pad(:, 2) ./ slices_interp) .* slices_interp ];
            % The upper limit is actually one slice too many, but it is
            % supposed to be zero anyway

            paddings = [ ...
                lims_blobs_w_orig(:, 1) - lims_blobs_w_pad(:, 1), ...
                lims_blobs_w_pad(:, 2) - lims_blobs_w_orig(:, 2) ];
        end

        function [lims_blobs_p_pad, lims_projs_p, paddings] = get_blob_p_lims(~, p_tab_lims, bl, slices_interp)
        % Finding limits for the shape functions