Gt6DReconstructionAlgorithmFactory.m

classdef Gt6DReconstructionAlgorithmFactory < handle
    properties
        data_type = 'single';

        vol_type = '3D';

        det_index = 1;

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

        use_predicted_scatter_ints = false;

        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, good_orients, blobs] = getReconstructionAlgo(self, sampler, num_interp, varargin)
            if (numel(sampler) == 1)
                [algo, good_orients, blobs] = self.getGrainReconstructionAlgo(sampler, num_interp, varargin{:});
            else
                [algo, good_orients, blobs] = self.getTwinReconstructionAlgo(sampler, num_interp, varargin{:});
            end
        end

        function [algo, good_orients, blobs] = getGrainReconstructionAlgo(self, sampler, num_interp, varargin)
            % Build Empty volumes
            ref_gr = sampler.get_reference_grain();
            proj = ref_gr.proj(self.det_index);

            spacing = mean([proj.vol_size_y, proj.vol_size_x, proj.vol_size_z]) * (self.rspace_oversize - 1);
            volume_size = round([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

            algo_params = self.get_algorithm_params(sampler, num_interp);

            good_orients = algo_params.good_orients;
            blobs = algo_params.blobs;

            % 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.
            if (~isempty(algo_params.geometries_ss))
                for ii_o = 1:numel(algo_params.offsets_ss)
                    offsets_ss = algo_params.offsets_ss{ii_o};
                    new_offsets_ss = offsets_ss{1};
                    cumulative_offs = numel(new_offsets_ss.proj);

                    for ii_ss = 2:numel(offsets_ss)
                        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);

                        num_projs = numel(offsets_ss{ii_ss}.proj);
                        for ii_p = 1:num_projs
                            offsets_ss{ii_ss}.proj(ii_p).sino_offset ...
                                = offsets_ss{ii_ss}.proj(ii_p).sino_offset + cumulative_offs;
                        end
                        new_offsets_ss.proj = cat(2, new_offsets_ss.proj, offsets_ss{ii_ss}.proj);

                        cumulative_offs = cumulative_offs + num_projs;
                    end
                    algo_params.offsets_ss{ii_o} = new_offsets_ss;
                end
            end

            algo = Gt6DBlobReconstructor(volume_size, blobs, ...
                'data_type', self.data_type, ...
                'geometries', algo_params.geometries, ...
                'offsets', algo_params.offsets, ...
                'ss_geometries', algo_params.geometries_ss, ...
                'ss_offsets', algo_params.offsets_ss, ...
                'use_astra_projectors', true, ...
                'volume_ss', self.volume_downscaling, ...
                'rspace_ss', self.rspace_super_sampling, ...
                'psf', algo_params.psf, ...
                varargin{:} );
        end

        function [algo, good_orients, blobs] = getTwinReconstructionAlgo(self, sampler, num_interp, varargin)
            num_grains = numel(sampler);

            % Build Empty volumes
            ref_gr = sampler(1).get_reference_grain();
            proj = ref_gr.proj(self.det_index);

            spacing = mean([proj.vol_size_y, proj.vol_size_x, proj.vol_size_z]) * (self.rspace_oversize - 1);
            volume_size = round([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

            % This vector will tell us where we will find in the blobs, the
            % ones that were selected among the included, because the
            % information we have in the twins about the spots that
            % coincide are about the included and not the selected
            included_pos_in_selected = zeros(size(proj.included));
            included_pos_in_selected(proj.selected) = 1:numel(find(proj.selected));

            blobs_to_not_shrink = cell(num_grains, 1);
            blobs_to_not_shrink{1} = false(size(included_pos_in_selected));
            for ii_g = 2:num_grains
                temp_gr = sampler(ii_g).get_reference_grain();
                shared = temp_gr.proj.shared_parent(sampler(ii_g).selected) > 0;
                shared_pos = temp_gr.proj.shared_parent(sampler(ii_g).selected);

                blobs_to_not_shrink{ii_g} = false(size(shared_pos));
                blobs_to_not_shrink{ii_g}(shared) = true;

                blobs_to_not_shrink{1}(included_pos_in_selected(shared_pos(shared))) = true;
            end

            for ii_g = 1:num_grains
                gr = sampler(ii_g).get_reference_grain();
                fprintf('%d) Grainid %d:\n', ii_g, gr.id)
                algo_params(ii_g) = self.get_algorithm_params(sampler(ii_g), num_interp, blobs_to_not_shrink{ii_g}); %#ok<AGROW>
                fprintf('\n')
            end

            blobs = cat(1, algo_params(:).blobs);

            psf = cat(1, algo_params(:).psf);

            good_orients = cat(1, algo_params(:).good_orients);

            geometries = cat(1, algo_params(:).geometries);
            geometries_ss = cat(1, algo_params(:).geometries_ss);

            offsets = algo_params(1).offsets;
            offsets_ss = algo_params(1).offsets_ss;
            base_offset_shift = numel(algo_params(1).blobs);
            for ii_g = 2:num_grains
                temp_gr = sampler(ii_g).get_reference_grain();
                shared = temp_gr.proj.shared_parent(sampler(ii_g).selected) > 0;
                shared_pos = temp_gr.proj.shared_parent(sampler(ii_g).selected);

                orient_offsets = algo_params(ii_g).offsets;
                for ii_o = 1:numel(orient_offsets)
                    to_be_summed = ~shared(orient_offsets{ii_o}.blob_offsets);

                    orient_offsets{ii_o}.blob_offsets(to_be_summed) ...
                        = orient_offsets{ii_o}.blob_offsets(to_be_summed) + base_offset_shift;

                    orient_offsets{ii_o}.blob_offsets(~to_be_summed) ...
                        = included_pos_in_selected(shared_pos(orient_offsets{ii_o}.blob_offsets(~to_be_summed)));

                    for ii_p = 1:numel(orient_offsets{ii_o}.proj)
                        ii_b = orient_offsets{ii_o}.proj(ii_p).blob_offset;
                        if (~shared(ii_b))
                            orient_offsets{ii_o}.proj(ii_p).blob_offset ...
                                = ii_b + base_offset_shift;
                        else
                            orient_offsets{ii_o}.proj(ii_p).blob_offset ...
                                = included_pos_in_selected(shared_pos(ii_b));
                        end
                    end
                end
                offsets = cat(1, offsets, orient_offsets);

                % Shifting projection to the concatenated blobs, unless that
                % blob was shared between parent and twin
                if (~isempty(geometries_ss))
                    orient_offsets_ss = algo_params(ii_g).offsets_ss;
                    for ii_o = 1:numel(orient_offsets_ss)
                        for ii_ss = 1:numel(orient_offsets_ss{ii_o})
                            to_be_summed = ~shared(orient_offsets{ii_o}{ii_ss}.blob_offsets);

                            orient_offsets_ss{ii_o}{ii_ss}.blob_offsets(to_be_summed) ...
                                = orient_offsets_ss{ii_o}{ii_ss}.blob_offsets(to_be_summed) + base_offset_shift;

                            orient_offsets_ss{ii_o}{ii_ss}.blob_offsets(~to_be_summed) ...
                                = included_pos_in_selected(shared_pos(orient_offsets_ss{ii_o}{ii_ss}.blob_offsets(~to_be_summed)));

                            for ii_p = 1:numel(orient_offsets{ii_o}{ii_ss}.proj)
                                ii_b = orient_offsets{ii_o}{ii_ss}.proj(ii_p).blob_offset;
                                if (~shared(ii_b))
                                    orient_offsets{ii_o}{ii_ss}.proj(ii_p).blob_offset ...
                                        = ii_b + base_offset_shift;
                                else
                                    orient_offsets{ii_o}{ii_ss}.proj(ii_p).blob_offset ...
                                        = included_pos_in_selected(shared_pos(ii_b));
                                end
                            end
                        end
                    end
                    offsets_ss = cat(1, offsets_ss, orient_offsets_ss);
                end

                base_offset_shift = base_offset_shift + numel(algo_params(ii_g).blobs);
            end

            algo = Gt6DBlobReconstructor(volume_size, blobs, ...
                'data_type', self.data_type, ...
                'geometries', geometries, ...
                'offsets', offsets, ...
                'ss_geometries', geometries_ss, ...
                'ss_offsets', offsets_ss, ...
                'use_astra_projectors', true, ...
                'volume_ss', self.volume_downscaling, ...
                'psf', psf, ...
                varargin{:} );
        end
    end

    methods (Access = protected)
        function algo_params = get_algorithm_params(self, sampler, num_interp, blobs_to_not_shrink)
            fprintf('Extracting blobs on detector: ')
            c = tic();
            sel_bls = sampler.selected;
            bl = sampler.bl(sel_bls);
            fprintf('\b\b (%f s).\n', toc(c));
            tot_blobs = numel(bl);
            sel_incl_indx = find(sel_bls);

            if (~exist('blobs_to_not_shrink', 'var') || isempty(blobs_to_not_shrink))
                blobs_to_not_shrink = false(size(bl));
            end

            % Should be decided automatically
            padding = 10;

            ref_gr = sampler.get_reference_grain();
            proj = ref_gr.proj(self.det_index);

            sel_allb_idxs = sampler.ondet(sampler.included(sel_bls));
            etas = ref_gr.allblobs.eta(sel_allb_idxs);
            omegas = ref_gr.allblobs.omega(sel_allb_idxs);

            slices_interp = self.compute_lorentz_slice_broadening(num_interp, etas);

            [orients, orients_ss] = sampler.get_orientations();
            orients = [orients{:}];
            % Array of structures that contain all the info relative to
            % each orientation, for each blob
            grains_props = self.get_selected_props(orients, sel_bls);

            % We should use sub-w interpolation, especially for num_inerp=1
            [w_tab, ref_ws] = self.get_w_tab(grains_props, omegas);

            [w_tab, slice_lims, lims, abs_lims, extreemes_blobs_w, extreemes_projs_w] ...
                = self.get_sub_blob_lims(w_tab, bl, ref_ws, slices_interp, padding);

            [sub_blob_slices, proj_coeffs] = self.get_sub_blobs(bl, slices_interp, padding);
            if (self.use_predicted_scatter_ints)
                sub_blob_slices = self.renormalize_blobs(sampler, sub_blob_slices);
            end

            [overflow_b, overflow_o] = self.compute_overflow(w_tab, extreemes_blobs_w);

            % We have to remove the orientations that don't project to any
            % blob
            tot_orientations = numel(grains_props);
            good_orients = overflow_o < tot_blobs;
            grains_props = grains_props(good_orients);
            w_tab = w_tab(:, good_orients);

            % We have to remove the blobs where no orientation projects
            % (keeping the total number of orientations intact for the
            % moment, because these numbers were computed with the previous
            % number of orientations).
            good_blobs = overflow_b < tot_orientations;
            if (numel(find(good_blobs)) < tot_blobs)
                fprintf(['WARNING! Some selected blobs were filtered out, ' ...
                    'due to no orientation projecting to the said blob:\n'])
                disp(find(~good_blobs)')
                fprintf('Position in the included vector:\n')
                incl_pos = false(size(sel_bls));
                incl_pos(sel_incl_indx(~good_blobs)) = true;
                disp(find(incl_pos)')

                titles = {'blobs img lims:', 'blobs img lims:', 'projs img lims:'};
                fprintf(' %15s %15s %15s\n', titles{:})
                fprintf('       %4d %4d       %4d %4d       %4d %4d\n', ...
                    [extreemes_blobs_w(~good_blobs, :), cat(1, bl(~good_blobs).bbwim), ...
                    round(extreemes_projs_w(~good_blobs, :))]');
            end
            bl = bl(good_blobs);
            sub_blob_slices = sub_blob_slices(good_blobs);
            proj_coeffs = proj_coeffs(good_blobs);
            etas = etas(good_blobs);
            omegas = omegas(good_blobs);
            slices_interp = slices_interp(good_blobs);
            slice_lims = slice_lims(good_blobs, :);
            abs_lims = abs_lims(good_blobs, :);
            lims = lims(good_blobs, :);
            w_tab = w_tab(good_blobs, :);
            extreemes_blobs_w = extreemes_blobs_w(good_blobs, :);
            extreemes_projs_w = extreemes_projs_w(good_blobs, :);
            overflow_b = overflow_b(good_blobs);

            tot_blobs = numel(find(good_blobs));
            sel_incl_indx = sel_incl_indx(good_blobs);

            if (~isempty(orients_ss))
                orients_props_ss = [orients_ss{:}];
                orients_props_ss = [orients_props_ss{:}];

                orients_props_ss = self.get_selected_props(orients_props_ss, sel_incl_indx);
                w_tab_ss = self.get_w_tab(orients_props_ss);
                [w_tab_ss, slice_lims_ss, lims_ss, abs_lims_ss, extreemes_blobs_w_ss, extreemes_projs_w_ss] ...
                    = self.get_sub_blob_lims(w_tab_ss, bl, ref_ws, slices_interp, padding);

                [overflow_ss_b, overflow_ss_o] = self.compute_overflow(w_tab_ss, extreemes_blobs_w_ss);

                good_ss_or = overflow_ss_o < tot_blobs;
                orients_props_ss = orients_props_ss(good_ss_or);
                w_tab_ss = w_tab_ss(:, good_ss_or);

                lims_tot = [ ...
                    min(lims(:, 1), lims_ss(:, 1)), ...
                    max(lims(:, 2), lims_ss(:, 2)) ];
                abs_lims_tot = [ ...
                    min(abs_lims(:, 1), abs_lims_ss(:, 1)), ...
                    max(abs_lims(:, 2), abs_lims_ss(:, 2)) ];
                slice_lims_tot = [ ...
                    min(slice_lims(:, 1), slice_lims_ss(:, 1)), ...
                    max(slice_lims(:, 2), slice_lims_ss(:, 2)) ];

                extreemes_blobs_w_tot = [ ...
                    min(extreemes_blobs_w(:, 1), extreemes_blobs_w_ss(:, 1)), ...
                    max(extreemes_blobs_w(:, 2), extreemes_blobs_w_ss(:, 2)) ];
                extreemes_projs_w_tot = [ ...
                    min(extreemes_projs_w(:, 1), extreemes_projs_w_ss(:, 1)), ...
                    max(extreemes_projs_w(:, 2), extreemes_projs_w_ss(:, 2)) ];
            else
                lims_tot = lims;
                abs_lims_tot = abs_lims;
                slice_lims_tot = slice_lims;

                overflow_ss_b = zeros(size(overflow_b));
                overflow_ss_o = zeros(size(overflow_o));

                extreemes_blobs_w_tot = extreemes_blobs_w;
                extreemes_projs_w_tot = extreemes_projs_w;
            end

            blobs = cell(size(sub_blob_slices));
            % Let's restrict blobs to the reached regions
            for ii = 1:tot_blobs
                fprintf('Blob: %03d (n. %03d) --> ', ii, sel_incl_indx(ii))
                if (overflow_b(ii) || overflow_ss_b(ii))
                    fprintf('Some orientations (%d, ss %d) will be missing.\n', ...
                        overflow_b(ii), overflow_ss_b(ii) )
                else
                    fprintf('Every orientation is projecting to the blob.\n')
                end

                if (slice_lims_tot(ii, 2) > size(sub_blob_slices{ii}, 3))
                    depth_blob = size(bl(ii).intm, 3);
                    fprintf('\n ERROR! numInterp: %d, lims: %d:%d, blob size: 1:%d, eta: %f, virtual slices 1:%d\n\n', ...
                        slices_interp(ii), lims_tot(ii, :), depth_blob, etas(ii), ...
                        depth_blob + mod(slices_interp(ii) - mod(depth_blob, slices_interp(ii)), ...
                        slices_interp(ii)) + 1)
                end

                fprintf(' - ExtremesBlobs %4d:%4d, ExtremesProjs %4d:%4d, Lims %2d:%2d, SliceLims %2d:%2d, Size 1:%2d ', ...
                    extreemes_blobs_w_tot(ii, :), ...
                    round(extreemes_projs_w_tot(ii, :)), ...
                    round(lims_tot(ii, :)), slice_lims_tot(ii, :), ...
                    size(sub_blob_slices{ii}, 3))

                blobs{ii} = sub_blob_slices{ii};
                if (~blobs_to_not_shrink(ii))
                    blobs{ii} = blobs{ii}(:, :, slice_lims_tot(ii, 1):slice_lims_tot(ii, 2));
                end
                blobs{ii} = permute(blobs{ii}, [1 3 2]);

                intensIncl = sum(sum(sum(blobs{ii}, 2)));
                intensExcl = sum(sum(sum(sub_blob_slices{ii}(:, :, [1:(slice_lims_tot(ii, 1)-1), (slice_lims_tot(ii, 2)+1):end]), 2)));
                fprintf('-> Intensity: included %3.3f, excluded %3.3f\n', ...
                    intensIncl, intensExcl);
            end

            fprintf('Computing blobs <-> sinograms coefficients: ')
            c = tic();
            [geometries, offsets] = self.makeSubBlobGeometries( ...
                grains_props, proj_coeffs, extreemes_blobs_w_tot, ...
                w_tab, abs_lims_tot, slice_lims_tot, blobs_to_not_shrink);
            fprintf('Done (%f s).\n', toc(c));

            num_geoms = numel(geometries);

            if (~isempty(orients_ss))
                geometries_ss = cell(num_geoms, 1);
                offsets_ss = cell(num_geoms, 1);
                fprintf('Computing blobs <-> supersampling_sinograms coefficients: ')
                orients_props_ss = orients_ss(good_orients);
                for ii = 1:num_geoms
                    num_chars = fprintf('%03d/%03d', ii, num_geoms);
                    orient_props_ss = orients_props_ss{ii};
                    orient_props_ss = [orient_props_ss{:}];
                    orient_props_ss = self.get_selected_props(orient_props_ss, sel_incl_indx);
                    w_tab_oss = self.get_w_tab(orient_props_ss);
                    w_tab_oss = self.fix_360_w_border(w_tab_oss, ref_ws);
                    [geometries_ss{ii}, offsets_ss{ii}] = self.makeSubBlobGeometries( ...
                        orient_props_ss, proj_coeffs, extreemes_blobs_w_tot, ...
                        w_tab_oss, abs_lims_tot, slice_lims_tot, blobs_to_not_shrink);
                    fprintf(repmat('\b', [1 num_chars]));
                end
                fprintf('Done (%f s).\n', toc(c));
            else
                geometries_ss = {};
                offsets_ss = {};
            end

            if (self.volume_downscaling > 1)
                for ii = 1:num_geoms
                    geometries{ii}(:, 4:12) = geometries{ii}(:, 4:12) / self.volume_downscaling;
                end
                if (~isempty(geometries_ss))
                    for ii = 1:num_geoms
                        geom_ss = geometries_ss{ii};
                        for ii_ss = 1:numel(geom_ss)
                            geom_ss{ii_ss}(:, 4:12) = geom_ss{ii_ss}(:, 4:12) / self.volume_downscaling;
                        end
                        geometries_ss{ii} = geom_ss;
                    end
                end
            end

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

            algo_params = struct( ...
                'blobs', {blobs}, ...
                'geometries', {geometries}, ...
                'offsets', {offsets}, ...
                'geometries_ss', {geometries_ss}, ...
                'offsets_ss', {offsets_ss}, ...
                'psf', {psf}, ...
                'good_orients', {good_orients'} );
        end

        function [sub_blob_slices, proj_coeffs] = get_sub_blobs(self, blobs, slices_interp, padding)
            num_blobs = numel(blobs);
            sub_blob_slices = cell(num_blobs, 1);
            sub_blob_coeffs = cell(num_blobs, 1);
            proj_coeffs = cell(num_blobs, 1);

            blob_sizes = zeros(num_blobs, 3);
            for ii = 1:num_blobs
                [blob_sizes(ii, 1), blob_sizes(ii, 2), blob_sizes(ii, 3)] = size(blobs(ii).intm);
            end

            % Number of interpolated slices should be at least enough
            % to accomodate all the slices of the blob
            nums_slices = ceil((blob_sizes(:, 3) -1) ./ slices_interp) +1 + 2 * padding;

            covered_slices_size = (nums_slices - 1) .* slices_interp + 1;

            fprintf('Creating sub-blob volumes: ')
            c = tic();
            for ii = 1:num_blobs
                num_chars = fprintf('%03d/%03d', ii, num_blobs);

                % Let's build the interpolated slices
                blob = blobs(ii).intm;
                blob(blob < 0) = 0;
                if (self.use_predicted_scatter_ints)
                    blob = blob .* blobs(ii).intensity;
                end

                sub_blob_slices{ii} = zeros(blob_sizes(ii, 1), blob_sizes(ii, 2), nums_slices(ii), self.data_type);
                for n = nums_slices(ii):-1:1
                    coeffs = abs( (slices_interp(ii) * (n-1) +1) - (1:covered_slices_size(ii)));
                    coeffs = 1 - coeffs / slices_interp(ii);

                    % these will be the blob slices that relate to the
                    % given ii interpolated slice
                    good_indxes = find(coeffs > 0);
                    % The corresponding coefficients
                    coeffs = coeffs(good_indxes);

                    blob_indexes = good_indxes - padding .* slices_interp(ii);
                    safe_indexes_pos = blob_indexes <= blob_sizes(ii, 3) & blob_indexes >= 1;
                    safe_indexes = blob_indexes(safe_indexes_pos);
                    safe_coeffs = coeffs(safe_indexes_pos);
                    for w = 1:numel(safe_coeffs)
                        sub_blob_slices{ii}(:, :, n) = ...
                            sub_blob_slices{ii}(:, :, n) + safe_coeffs(w) .* blob(:, :, safe_indexes(w));
                    end
                    % Save the metadata
                    sub_blob_coeffs{ii}(n) = struct('indexes', good_indxes, 'coeffs', coeffs);
                end

                % Perpare next metadata structure
                proj_coeffs{ii} = struct( ...
                    'slices', arrayfun(@(x){[]}, 1:covered_slices_size(ii)), ...
                    'coeffs', arrayfun(@(x){[]}, 1:covered_slices_size(ii)) );

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

            fprintf('\b\b (%f s), transforming coeffs..', toc(c));
            c = tic();

            % Let's transform the coeffs to something that is easy to
            % assign for projections:
            % Per each blob we want to easily find out what the projections
            % project to. The metadata from sub_blob_coeffs is the other
            % way around.
            for ii = 1:num_blobs
                num_slices = numel(sub_blob_coeffs{ii});

                for n = num_slices:-1:1
                    indexes = sub_blob_coeffs{ii}(n).indexes;
                    coeffs = sub_blob_coeffs{ii}(n).coeffs;
                    for m = numel(indexes):-1:1
                        proj_coeffs{ii}(indexes(m)).slices = ...
                            [ proj_coeffs{ii}(indexes(m)).slices, n ];
                        proj_coeffs{ii}(indexes(m)).coeffs = ...
                            [ proj_coeffs{ii}(indexes(m)).coeffs, coeffs(m) ];
                    end
                end
            end
            fprintf('\b\b (%f s).\n', toc(c));
        end

        function [w_tab, ref_w] = get_w_tab(self, orient_props, ref_ws)
            if (~isfield(self.parameters.labgeo, 'omstep'))
                omega_step = gtGetOmegaStepDeg(self.parameters);
            else
                omega_step = self.parameters.labgeo.omstep;
            end
            % Tabular of omegas: horizontal direction the n geometries,
            % vertical direction the omegas relative to each blob
            real_w_tab = [orient_props(:).omega];
            % Images numbers
            w_tab = real_w_tab / omega_step;

            if (nargout > 1)
                if (exist('ref_ws', 'var'))
                    ref_w = ref_ws / omega_step;
                else
                    ref_w = zeros(size(w_tab, 1), 1);
                end
            end
        end

        function w_tab = fix_360_w_border(self, w_tab, ref_int_ws)
            num_images = gtGetTotNumberOfImages(self.parameters);

            % Let's treat those blobs at the w edge 360->0
            % (from the sampled orientations perspective)
            opposite_int_ws = mod(ref_int_ws + num_images / 2, num_images);
            gt_ref_int = opposite_int_ws > num_images / 2;
            lt_ref_int = ~gt_ref_int;
            opposite_int_ws_repmat = opposite_int_ws(:, ones(1, size(w_tab, 2)));

            w_tab(gt_ref_int, :) = w_tab(gt_ref_int, :) ...
                - num_images .* (w_tab(gt_ref_int, :) > opposite_int_ws_repmat(gt_ref_int, :));
            w_tab(lt_ref_int, :) = w_tab(lt_ref_int, :) ...
                + num_images .* (w_tab(lt_ref_int, :) < opposite_int_ws_repmat(lt_ref_int, :));
        end

        function [w_tab, sliceLims, lims, absLims, extreemes_blobs_w, extreemes_projs_w] ...
                = get_sub_blob_lims(self, w_tab, bl, ref_int_ws, slices_interp, padding)
            w_tab = self.fix_360_w_border(w_tab, ref_int_ws);

            % Let's add extreme values
            extreemes_projs_w = [ ...
                min(w_tab, [], 2), ...
                max(w_tab, [], 2) ];

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

            blobs_sizes = extreemes_blobs_w(:, 2) - extreemes_blobs_w(:, 1) +1;

            extreemes_blobs_w(:, 2) = extreemes_blobs_w(:, 1) + (blobs_sizes-1) ...
                + mod(slices_interp - mod((blobs_sizes-1), slices_interp), slices_interp);

            % Adding the padding
            extreemes_blobs_w = [ ...
                extreemes_blobs_w(:, 1) - padding .* slices_interp, ...
                extreemes_blobs_w(:, 2) + padding .* slices_interp ];

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

            % Blob limits in the full omega stack
            absLims = [ ...
                max([extreemes_projs_w(:, 1), extreemes_blobs_w(:, 1)], [], 2), ...
                min([extreemes_projs_w(:, 2), extreemes_blobs_w(:, 2)], [], 2) ];
            % blob limits in its volume size
            lims = absLims - extreemes_blobs_w(:, [1 1]) + 1;

            sliceLims = (lims -1) ./ slices_interp(:, [1 1]);
            sliceLims = [floor(sliceLims(:, 1)), ceil(sliceLims(:, 2))] + 1;
        end

        function volume_size = get_volume_size(self, gr_center, Ws, bls)
            % At this point we kinda lost the original blob size
            % information! so we get it from the mask ;)
            BBs(:, [1, 3]) = cat(1, bls(:).bbuim);
            BBs(:, [2, 4]) = cat(1, bls(:).bbvim);

            BBs = [BBs(:, 1:2), (BBs(:, 3:4) - BBs(:, 1:2) + 1)];

            verts = gtFwdSimComputeCircumscribingPolyhedron(gr_center, Ws, BBs, self.parameters, self.det_index);

            volume_size = round(max(verts, [], 1) - min(verts, [], 1) / self.parameters.fsim.oversize);
        end

        function [geometries, offsets] = makeSubBlobGeometries( self, ...
                grains_props, proj_coeffs, extreemes_blobs_w, w_tab, abs_lims, slice_lims, blobs_to_not_shrink)
            % ASTRA Geometry
            if (isfield(grains_props, 'detector'))
                geometries = arrayfun(@(x){x.detector(self.det_index).proj_geom}, grains_props);
            else
                geometries = {grains_props(:).proj_geom};
            end
            geometries = reshape(geometries, [], 1);

            num_geoms = numel(geometries);

            w_int_tab(:, :, 1) = floor(w_tab);
            w_int_tab(:, :, 2) = ceil(w_tab);
            w_int_coeffs(:, :, 1) = 1 - (w_int_tab(:, :, 2) - w_tab);
            w_int_coeffs(:, :, 2) = 1 - w_int_coeffs(:, :, 1);

            slice_pos_in_blob = w_int_tab - repmat(extreemes_blobs_w(:, 1), [1 num_geoms 2]) + 1;

            % Assign offsets
            offsets = cell(num_geoms, 1);
            for ii = 1:num_geoms
                valid_proj_dirs = ...
                    (w_tab(:, ii) >= abs_lims(:, 1)) ...
                    & (w_tab(:, ii) <= abs_lims(:, 2));

                if (all(~valid_proj_dirs))
                    warning('Geometry:wrong_orientation', ...
                        'This orientation is useless! %d', ii)
                end

                valid_proj_indxes = find(valid_proj_dirs)';
                num_valid_projs = numel(valid_proj_indxes);

                offsets{ii}.('sino_offsets') = [];
                offsets{ii}.('blob_offsets') = [];
                offsets{ii}.('proj_offsets') = [];
                offsets{ii}.('proj_coeffs')  = [];

                offsets{ii}.proj(num_valid_projs) = struct( ...
                    'sino_offset', [], 'blob_offset', [], ...
                    'proj_offsets', [], 'proj_coeffs', [] );
                for jj = 1:num_valid_projs
                    n = valid_proj_indxes(jj);

                    slices_1 = [proj_coeffs{n}(slice_pos_in_blob(n, ii, 1)).slices];
                    coeffs_1 = [proj_coeffs{n}(slice_pos_in_blob(n, ii, 1)).coeffs] * w_int_coeffs(n, ii, 1);

                    slices_2 = [proj_coeffs{n}(slice_pos_in_blob(n, ii, 2)).slices];
                    coeffs_2 = [proj_coeffs{n}(slice_pos_in_blob(n, ii, 2)).coeffs] * w_int_coeffs(n, ii, 2);

                    % ugly but it avoids a lot of compuation later
                    if (numel(slices_1) == 2 && numel(slices_2) == 2 && all(slices_1 == slices_2))
                        slices = slices_1;
                        coeffs = coeffs_1 + coeffs_2;
                    elseif (numel(slices_1) == 1 && numel(slices_2) == 2)
                        ind = find(slices_2 == slices_1);
                        slices = slices_2;
                        coeffs = coeffs_2;
                        coeffs(ind) = coeffs(ind) + coeffs_1;
                    elseif (numel(slices_1) == 2 && numel(slices_2) == 1)
                        ind = find(slices_2 == slices_1);
                        slices = slices_1;
                        coeffs = coeffs_1;
                        coeffs(ind) = coeffs(ind) + coeffs_2;
                    else
                        slices = [slices_1, slices_2];
                        coeffs = [coeffs_1, coeffs_2];
                    end

                    % Correcting w-slices indices for blob truncation
                    if (~blobs_to_not_shrink(n))
                        slices = slices - slice_lims(n, 1) +1;
                    end

                    offsets{ii}.proj(jj) = struct( ...
                        'sino_offset', jj, 'blob_offset', n, ...
                        'proj_offsets', slices, 'proj_coeffs', coeffs );

                    offsets{ii}.('sino_offsets') = [offsets{ii}.('sino_offsets'), jj(1, ones(numel(slices), 1))];
                    offsets{ii}.('blob_offsets') = [offsets{ii}.('blob_offsets'), n(1, ones(numel(slices), 1))];
                    offsets{ii}.('proj_offsets') = [offsets{ii}.('proj_offsets'), slices];
                    offsets{ii}.('proj_coeffs')  = [offsets{ii}.('proj_coeffs'),  coeffs];
                end

                geometries{ii} = double(geometries{ii}(valid_proj_dirs, :));
            end
        end

        function blobs = renormalize_blobs(self, sampler, blobs)
            ref_gr = sampler.get_reference_grain();
            if (isfield(self.parameters.cryst(ref_gr.phaseid), 'int_exp'))
                cryst_scatter_ints = self.parameters.cryst(ref_gr.phaseid).int_exp;
            else
                cryst_scatter_ints = self.parameters.cryst(ref_gr.phaseid).int;
            end

            sel = sampler.ondet(sampler.included(sampler.selected));
            thetatypes = ref_gr.allblobs.thetatype(sel);
            scatter_ints = cryst_scatter_ints(thetatypes);
            lorentz_fact = ref_gr.allblobs.lorentzfac(sel);
            for ii_b = 1:numel(blobs)
                blobs{ii_b} = blobs{ii_b} ./ scatter_ints(ii_b) ./ lorentz_fact(ii_b);
            end
        end
    end

    methods (Access = protected, Static)
        function grains_props = get_selected_props(grains, blobs_on_det)
            for ii = numel(grains):-1:1
                grains(ii) = gtGrainAllblobsFilterOrder(grains(ii), ...
                    [], blobs_on_det);
            end
            grains_props = [grains(:).allblobs];
        end

        function allblobs = selectBlobsFromGrain(allblobs, blobs_on_det)
            names = fieldnames(allblobs);
            for ii = 1:numel(names)
                name = names{ii};
                if (strcmpi(name, 'srot'))
                    allblobs.(name) = allblobs.(name)(:, :, blobs_on_det);
                else
                    allblobs.(name) = allblobs.(name)(blobs_on_det, :);
                end
            end
        end

        function slices_interp = compute_lorentz_slice_broadening(num_interp, etas)
            lorentz_factor = 1 ./ abs(sind(etas));
            % Let's try to not penalize rounding errors and being a bit
            % tollerant with 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
            slices_interp = ceil(num_interp .* (lorentz_factor - 0.25));
        end

        function [overflow_b, overflow_o] = compute_overflow(w_int_tab, extreemes_blobs_w)
            num_orient = size(w_int_tab, 2);
            ones_orient = ones(1, num_orient);
            lower_lims = w_int_tab < extreemes_blobs_w(:, 1 .* ones_orient);
            upper_lims = w_int_tab > extreemes_blobs_w(:, 2 .* ones_orient);
            overflow_b = sum(lower_lims, 2) + sum(upper_lims, 2);
            overflow_o = sum(lower_lims, 1) + sum(upper_lims, 1);
        end
    end
end