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

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

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

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

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

            % computing the needed extension for the interpolation
            lims_blobs_p_pad = [ ...
                floor(lims_blobs_p_pad(:, 1) ./ slices_interp) .* slices_interp, ...
                ceil(lims_blobs_p_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_p_pad(:, 1), ...
                lims_blobs_p_pad(:, 2) - lims_blobs_w_orig(:, 2) ];
        end

        function geoms = compute_proj_geom_numinterp(self, projs_uv_lims, ref_gr, ors_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));

            num_ors = numel(ors_allblobs);
            om = zeros(numel(ab_sel), num_ors);

            for ii_o = 1:num_ors
                om(:, ii_o) = ors_allblobs(ii_o).omega(ab_sel);
            end
            if (isfield(self.parameters.acq, 'correct_sample_shifts'))
                [~, shifts_sam] = gtMatchGetSampleShifts(self.parameters, om(:));
            else
                shifts_sam = zeros(numel(om(:)), 3);
            end
            geoms = cell(num_ors, 1);
            for ii_o = 1:num_ors
                dvec_sam = ors_allblobs(ii_o).dvecsam(ab_sel, :);
                rot_w_l2s = ors_allblobs(ii_o).srot(:, :, ab_sel);

                % Dimensions are: <blobs x uv x [min, max]>
                proj_or_lims = projs_uv_lims(:, :, :, ii_o);
                bb_bls_uv = [proj_or_lims(:, :, 1), proj_or_lims(:, :, 2) - proj_or_lims(:, :, 1) + 1];
                shifts_sam_ii_o = shifts_sam((ii_o - 1) * numel(ab_sel) +1 : ii_o * numel(ab_sel), :);
                proj_geom = gtGeoProjForReconstruction(...
                    dvec_sam, rot_w_l2s, ref_gr.center, bb_bls_uv, [], ...
                    self.parameters.detgeo(det_ind), ...
                    self.parameters.labgeo, ...
                    self.parameters.samgeo, ...
                    self.parameters.recgeo(self.det_index(1)), ...
                    'ASTRA_grain', shifts_sam_ii_o);

                if (self.volume_downscaling > 1)
                    proj_geom(:, 4:12) = proj_geom(:, 4:12) / self.volume_downscaling;
                end

                geoms{ii_o} = double(proj_geom);
            end
        end

        function offsets = compute_proj_coeffs_numinterp(~, lims_blobs_w, w_tab, slices_interp, shifts_uv)
            % Dimensions of 'shifts_uv' are: <blobs x uv x orientations>

            num_blobs = size(w_tab, 1);
            num_orients = size(w_tab, 2);

            rescaled_w_tab = w_tab ./ slices_interp(:, ones(num_orients, 1));

            w_int_tab(:, :, 1) = floor(rescaled_w_tab);
            w_int_tab(:, :, 2) = ceil(rescaled_w_tab);

            slice_coeffs_in_blob(:, :, 1) = w_int_tab(:, :, 2) - rescaled_w_tab;
            slice_coeffs_in_blob(:, :, 2) = 1 - slice_coeffs_in_blob(:, :, 1);

            rescaled_low_lims_blobs_w = lims_blobs_w(:, 1) ./ slices_interp;
            slice_pos_in_blob = w_int_tab - repmat(rescaled_low_lims_blobs_w, [1 num_orients 2]) + 1;

            slice_coeffs_in_blob = permute(slice_coeffs_in_blob, [3 1 2]);
            slice_coeffs_in_blob = reshape(slice_coeffs_in_blob, 2 * num_blobs, num_orients);
            slice_pos_in_blob = permute(slice_pos_in_blob, [3 1 2]);
            slice_pos_in_blob = reshape(slice_pos_in_blob, 2 * num_blobs, num_orients);

            shifts_uv = double(shifts_uv);
            slice_coeffs_in_blob = double(slice_coeffs_in_blob);
            slice_pos_in_blob = double(slice_pos_in_blob);

            inds = 1:num_blobs;
            inds = reshape(inds([1 1], :), 1, []);

            % Assign offsets
            offsets = cell(num_orients, 1);
            for ii_o = 1:num_orients
                offsets{ii_o} = struct( ...
                    'sino_offsets', inds, ...
                    'blob_offsets', inds, ...
                    'proj_offsets', slice_pos_in_blob(:, ii_o)', ...
                    'proj_coeffs', slice_coeffs_in_blob(:, ii_o)', ...
                    'proj_shifts_uv', shifts_uv(inds, :, ii_o)');
            end
        end

        function [geometries, offsets] = compute_proj_geom_shapefuncs_w(...
                self, projs_uv_lims, ref_gr, ors_allblobs, lims_blobs_w, sf_w, blobs_w_interp, shifts_uv, det_ind)
        % We here produce ASTRA's projection geometry and blob assempling
        % coefficients
            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_blobs = numel(sf_w{1});
            num_orients = numel(ors_allblobs);

            rot_comp_w = gtMathsRotationMatrixComp(self.parameters.labgeo.rotdir', 'col');

            om_step = gtAcqGetOmegaStep(self.parameters, det_ind);

            geometries = cell(num_orients, 1);
            offsets = cell(num_orients, 1);
            for ii_o = 1:num_orients
                num_chars = fprintf('%03d/%03d', ii_o, num_orients);

                sf = sf_w{ii_o};

                % Assign offsets
                ind_of_blob = cell(1, num_blobs);
                slice_coeff = cell(1, num_blobs);
                ws = cell(1, num_blobs);
                for ii_b = 1:num_blobs
                    cs = reshape(sf(ii_b).intm, 1, []);
                    iws = sf(ii_b).bbwim(1):blobs_w_interp(ii_b):sf(ii_b).bbwim(2);

                    valid_coeffs = cs > eps('single');
                    slice_coeff{ii_b} = cs(valid_coeffs);
                    ws{ii_b} = iws(valid_coeffs);

                    ind_of_blob{ii_b} = ii_b * ones(1, sum(valid_coeffs));
                end

                ws = [ws{:}];
                pos_in_sino = 1:numel(ws);
                ind_of_blob = double([ind_of_blob{:}]);
                pos_in_blob = double((ws - lims_blobs_w(ind_of_blob, 1)') ./ blobs_w_interp(ind_of_blob, 1)' + 1);
                slice_coeff = double([slice_coeff{:}]);

                offsets{ii_o} = struct( ...
                    'sino_offsets', pos_in_sino, ...
                    'blob_offsets', ind_of_blob, ...
                    'proj_offsets', pos_in_blob, ...
                    'proj_coeffs', slice_coeff, ...
                    'proj_shifts_uv', shifts_uv(ind_of_blob, :, ii_o)' );

                % Computing geometries, but from dvecsam that have been
                % computed from scratch for each W associated with each W
                % of a shape-function
                rot_w_l2s = gtMathsRotationTensor(ws' .* om_step, rot_comp_w);

                dvec_lab = ors_allblobs(ii_o).dvec(ab_sel(ind_of_blob), :);
                dvec_sam = gtGeoLab2Sam(dvec_lab, rot_w_l2s, ...
                    self.parameters.labgeo, self.parameters.samgeo, true);

                % accounting for sample shifts relative to rotation axis
                if isfield(self.parameters.acq, 'correct_sample_shifts')
                    [~, shifts_sam] = gtMatchGetSampleShifts(self.parameters, ors_allblobs(ii_o).omega(ab_sel(ind_of_blob)));
                else
                    shifts_sam = zeros(numel(ors_allblobs(ii_o).omega(ab_sel(ind_of_blob))), 3);
                end

                % Dimensions are: <blobs x uv x [min, max]>
                proj_or_lims = projs_uv_lims(:, :, :, ii_o);
                bb_bls_uv = [proj_or_lims(:, :, 1), proj_or_lims(:, :, 2) - proj_or_lims(:, :, 1) + 1];

                proj_geom = gtGeoProjForReconstruction(...
                    dvec_sam, rot_w_l2s, ref_gr.center, ...
                    bb_bls_uv(ind_of_blob, :), [], ...
                    self.parameters.detgeo(det_ind), ...
                    self.parameters.labgeo, ...
                    self.parameters.samgeo, ...
                    self.parameters.recgeo(self.det_index(1)), ...
                    'ASTRA_grain', shifts_sam);

                if (self.volume_downscaling > 1)
                    proj_geom(:, 4:12) = proj_geom(:, 4:12) / self.volume_downscaling;
                end

                geometries{ii_o} = double(proj_geom);

                fprintf(repmat('\b', [1 num_chars]));
            end
        end

        function [geometries, offsets] = compute_proj_geom_shapefuncs_nw(...
                self, projs_uv_lims, ref_gr, ors_allblobs, lims_blobs_w, sf_nw, blobs_w_interp, shifts_uv, det_ind)
        % We here produce ASTRA's projection geometry and blob assempling
        % coefficients
            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_blobs = numel(sf_nw{1});
            num_orients = numel(ors_allblobs);

            rot_comp_w = gtMathsRotationMatrixComp(self.parameters.labgeo.rotdir', 'col');

            om_step = gtAcqGetOmegaStep(self.parameters, det_ind);

            geometries = cell(num_orients, 1);
            offsets = cell(num_orients, 1);
            for ii_o = 1:num_orients
                num_chars = fprintf('%03d/%03d', ii_o, num_orients);

                sf = sf_nw{ii_o};

                bl_bbsizes = cat(1, sf.bbsize);

                % Assign offsets
                ind_of_blob = cell(1, num_blobs);
                slice_coeff = cell(1, num_blobs);
                ns = cell(1, num_blobs);
                ws = cell(1, num_blobs);
                for ii_b = 1:num_blobs
                    cs = sf(ii_b).intm;
                    cs = reshape(cs, 1, []);
                    % Omegas of each column of the shape-function
                    iws = sf(ii_b).bbwim(1):blobs_w_interp(ii_b):sf(ii_b).bbwim(2);

                    % Etas of each row of the shape-function
                    ins = linspace(sf(ii_b).bbnim(1), sf(ii_b).bbnim(2), bl_bbsizes(ii_b, 1));

                    ins = reshape(ins, [], 1);
                    ins = ins(:, ones(numel(iws), 1));
                    ins = reshape(ins, [], 1);

                    iws = iws(ones(sf(ii_b).bbsize(1), 1), :);
                    iws = reshape(iws, [], 1);

                    % ins and iws have been expanded to be of the same
                    % size, but there is some redundancy going on there.
                    % This could impact on the perforance of the projection
                    % geometry determination later in this same function,
                    % but this is the only way to keep consistency with the
                    % sinogram<->blobs transformation coefficients.
                    valid_coeffs = cs > eps('single');
                    ns{ii_b} = ins(valid_coeffs);
                    ws{ii_b} = iws(valid_coeffs);
                    slice_coeff{ii_b} = cs(valid_coeffs);

                    ind_of_blob{ii_b} = ii_b * ones(1, sum(valid_coeffs));
                end

                ws = cat(1, ws{:});
                ns = cat(1, ns{:});
                pos_in_sino = 1:numel(ws);
                ind_of_blob = double([ind_of_blob{:}]);
                pos_in_blob = double(((ws - lims_blobs_w(ind_of_blob, 1)) ./ blobs_w_interp(ind_of_blob, 1))' + 1);
                slice_coeff = double([slice_coeff{:}]);

                offsets{ii_o} = struct( ...
                    'sino_offsets', pos_in_sino, ...
                    'blob_offsets', ind_of_blob, ...
                    'proj_offsets', pos_in_blob, ...
                    'proj_coeffs', slice_coeff, ...
                    'proj_shifts_uv', shifts_uv(ind_of_blob, :, ii_o)' );

                % Computing geometries
                rot_w_l2s = gtMathsRotationTensor(ws .* om_step, rot_comp_w);

                % Building the projection geometry from the etas and omegas
                t = ors_allblobs(ii_o).theta(ab_sel(ind_of_blob));
                pl_lab = gtGeoPlLabFromThetaEta(t, ns, self.parameters.labgeo);
                dvec_lab = gtFedPredictDiffVecMultiple(pl_lab', self.parameters.labgeo.beamdir')';
                dvec_sam = gtGeoLab2Sam(dvec_lab, rot_w_l2s, ...
                    self.parameters.labgeo, self.parameters.samgeo, true);

                % accounting for sample shifts relative to rotation axis
                if isfield(self.parameters.acq, 'correct_sample_shifts')
                    [~, shifts_sam] = gtMatchGetSampleShifts(self.parameters, ors_allblobs(ii_o).omega(ab_sel(ind_of_blob)));
                else
                    shifts_sam = zeros(numel(ors_allblobs(ii_o).omega(ab_sel(ind_of_blob))), 3);
                end

                % Dimensions are: <blobs x uv x [min, max]>
                proj_or_lims = projs_uv_lims(:, :, :, ii_o);
                bb_bls_uv = [proj_or_lims(:, :, 1), proj_or_lims(:, :, 2) - proj_or_lims(:, :, 1) + 1];

                proj_geom = gtGeoProjForReconstruction(...
                    dvec_sam, rot_w_l2s, ref_gr.center, ...
                    bb_bls_uv(ind_of_blob, :), [], ...
                    self.parameters.detgeo(det_ind), ...
                    self.parameters.labgeo, ...
                    self.parameters.samgeo, ...
                    self.parameters.recgeo(self.det_index(1)), ...
                    'ASTRA_grain', shifts_sam);

                if (self.volume_downscaling > 1)
                    proj_geom(:, 4:12) = proj_geom(:, 4:12) / self.volume_downscaling;
                end

                geometries{ii_o} = double(proj_geom);

                fprintf(repmat('\b', [1 num_chars]));
            end
        end

        function [scatter_ints_avg, scatter_ints_ors] = get_scattering_intensities(self, ref_or, ors_allbls, det_ind)
            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 (isfield(self.parameters.cryst(ref_or.phaseid), 'int_exp'))
                fprintf(' - Using experimental intensities for each family (from: parameters.cryst.int_exp)\n')
                cryst_scatter_ints = self.parameters.cryst(ref_or.phaseid).int_exp;
            else
                fprintf(' - Using theoretically predicted intensities for each family\n')
                cryst_scatter_ints = self.parameters.cryst(ref_or.phaseid).int;
            end
            cryst_scatter_ints = cryst_scatter_ints ./ mean(cryst_scatter_ints);

            thetatypes = ref_or.allblobs(det_ind).thetatype(ab_sel);
            scatter_ints = cryst_scatter_ints(thetatypes);
            scatter_ints = reshape(scatter_ints, [], 1);

            if (self.use_predicted_scatter_ints > 1)
                try
                    fprintf(' - Applying incoming beam intensity corrections\n')
                    [beam_intensities, refs] = gtGrainComputeIncomingBeamIntensity(ref_or, self.parameters, det_ind);
                    beam_in_ints_avg = beam_intensities(sel);
                    beam_in_ints_avg = beam_in_ints_avg / mean(beam_in_ints_avg);
                catch mexc
                    gtPrintException(mexc, 'Reverting to non-corrected beam intensities')
                    beam_in_ints_avg = 1;
                end

                try
                    fprintf(' - Applying beam attenuation corrections\n')
                    [atts_tot, ~, abs_vol] = gtGrainComputeBeamAttenuation(ref_or, self.parameters, det_ind);

                    beam_out_ints_avg = atts_tot(sel);
                catch mexc
                    gtPrintException(mexc, 'Reverting to non-corrected beam attenuations')
                    beam_out_ints_avg = 1;
                end
            else
                beam_in_ints_avg = 1;
                beam_out_ints_avg = 1;
            end

            scatter_ints_avg = scatter_ints .* ref_or.allblobs(det_ind).lorentzfac(ab_sel) .* beam_in_ints_avg .* beam_out_ints_avg;

            fprintf(' - Computing each orientation''s corrections..')
            c = tic();
            if (self.use_predicted_scatter_ints > 2 ...
                    && (numel(beam_in_ints_avg) > 1 || numel(beam_out_ints_avg) > 1))
                num_ors = numel(ors_allbls);
                scatter_ints_ors = scatter_ints(:, ones(num_ors, 1));
                fprintf('\b\b: ')
                for ii_o = 1:num_ors
                    if (mod(ii_o, 10) == 1)
                        num_chars = fprintf('%05d/%05d', ii_o, num_ors);
                    end
                    or = ref_or;
                    or.allblobs = ors_allbls(ii_o);
                    if (numel(beam_in_ints_avg) > 1)
                        try
                            [beam_in_int_or, refs] = gtGrainComputeIncomingBeamIntensity(or, self.parameters, det_ind, refs);
                            beam_in_int_or = beam_in_int_or(sel);
                            beam_in_int_or = beam_in_int_or / mean(beam_in_int_or);
                        catch mexc
                            gtPrintException(mexc, 'Reverting to non-corrected beam intensities')
                            beam_in_int_or = 1;
                        end
                    else
                        beam_in_int_or = 1;
                    end

                    if (numel(beam_out_ints_avg) > 1)
                        try
                            atts_tot = gtGrainComputeBeamAttenuation(or, self.parameters, det_ind, abs_vol);

                            beam_out_int_or = atts_tot(sel);
                        catch mexc
                            gtPrintException(mexc, 'Reverting to non-corrected beam intensities')
                            beam_out_int_or = 1;
                        end
                    else
                        beam_out_int_or = 1;
                    end

                    scatter_ints_ors(:, ii_o) = scatter_ints_ors(:, ii_o) .* ors_allbls(ii_o).lorentzfac(ab_sel, :) .* beam_in_int_or .* beam_out_int_or;

                    if (mod(ii_o, 10) == 1)
                        fprintf(repmat('\b', [1 num_chars]));
                    end
                end
            else
                lfacs = cat(2, ors_allbls(:).lorentzfac);
                scatter_ints_ors = scatter_ints .* beam_in_ints_avg .* beam_out_ints_avg;
                scatter_ints_ors = bsxfun(@times, scatter_ints_ors, lfacs(ab_sel, :));
            end
            fprintf('\b\b: Done in %g seconds.\n', toc(c))
        end

        function or_groups = get_orientation_groups(~, sampler)
            num_regions = numel(sampler);
            or_groups = zeros(num_regions, 2);
            or_groups(1, 1) = 1;
            for ii_o_reg = 2:num_regions
                or_groups(ii_o_reg, 1) = or_groups(ii_o_reg-1) + sampler(ii_o_reg-1).get_total_num_orientations();
            end
            or_groups(1:end-1, 2) = or_groups(2:end, 1) - 1;
            or_groups(end, 2) = or_groups(end, 1) + sampler(end).get_total_num_orientations() - 1;
        end
    end

    methods (Access = protected, Static)
        function new_offsets_ss = merge_offsets(offsets_ss)
            new_offsets_ss = offsets_ss{1};
            cumulative_offs = max(offsets_ss{1}.sino_offsets);

            num_sub_orients = numel(offsets_ss);

            for ii_ss = 2:num_sub_orients
                new_offsets_ss.sino_offsets ...
                    = cat(2, new_offsets_ss.sino_offsets, offsets_ss{ii_ss}.sino_offsets + cumulative_offs);
                new_offsets_ss.blob_offsets ...
                    = cat(2, new_offsets_ss.blob_offsets, offsets_ss{ii_ss}.blob_offsets);
                new_offsets_ss.proj_offsets ...
                    = cat(2, new_offsets_ss.proj_offsets, offsets_ss{ii_ss}.proj_offsets);
                new_offsets_ss.proj_coeffs ...
                    = cat(2, new_offsets_ss.proj_coeffs, offsets_ss{ii_ss}.proj_coeffs);
                new_offsets_ss.proj_shifts_uv ...
                    = cat(2, new_offsets_ss.proj_shifts_uv, offsets_ss{ii_ss}.proj_shifts_uv);

                num_projs = max(offsets_ss{ii_ss}.sino_offsets);
                cumulative_offs = cumulative_offs + num_projs;
            end

            new_offsets_ss.proj_coeffs = new_offsets_ss.proj_coeffs / num_sub_orients;
        end
    end
end
