classdef Gt6DReconstructionAlgorithmFactory < handle
    properties
        data_type = 'single';
    end

    methods (Access = public)
        function [algo, blobs] = getBlobReconstructionAlgo(self, bl, selected_bls, sampler, vol_type)
            fprintf('Extracting blobs on detector: ')
            c = tic();
            bl = bl(selected_bls);
            fprintf('\b\b (%f s).\n', toc(c));
            tot_blobs = numel(bl);

            [grains, grains_ss] = sampler.get_orientations();
            grains = [grains{:}];
            % Array of structures that contain all the info relative to
            % each orientation, for each blob
            grains_props = Gt6DReconstructionAlgorithmFactory.get_selected_props(grains, selected_bls);
            % Tabular of omegas: horizontal direction the n geometries,
            % vertical direction the omegas relative to each blob
            omega_int_tab = Gt6DReconstructionAlgorithmFactory.getOmegaIntegerTab(grains_props, sampler.parameters);
%             [omega_int_tab, omega_int_coeffs] = Gt6DReconstructionAlgorithmFactory.getOmegaIntegersAndCoeffs(grainsProps, sampler.parameters);

            [lims, abs_lims, ~, extremesOfBlobs, extremesOmIntTab] ...
                = Gt6DReconstructionAlgorithmFactory.getBlobLims(omega_int_tab, bl);

            [overflow_b, overflow_o] = Gt6DReconstructionAlgorithmFactory.compute_overflow(omega_int_tab, extremesOfBlobs);

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

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

                grains_props_ss = Gt6DReconstructionAlgorithmFactory.get_selected_props(grains_ss, selected_bls);
                omega_int_tab_ss = Gt6DReconstructionAlgorithmFactory.getOmegaIntegerTab(grains_props_ss, sampler.parameters);
                [lims_ss, abs_lims_ss, ~, extremesOfBlobs_ss, extremesOmIntTab_ss] ...
                    = Gt6DReconstructionAlgorithmFactory.getBlobLims(omega_int_tab_ss, bl);

                [overflow_ss_b, overflow_ss_o] = Gt6DReconstructionAlgorithmFactory.compute_overflow(omega_int_tab_ss, extremesOfBlobs_ss);

                good_ss_or = overflow_ss_o < tot_blobs;
                grains_props_ss = grains_props_ss(good_ss_or);
                omega_int_tab_ss = omega_int_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))];

                extremesOfBlobs_tot = [min(extremesOfBlobs(:, 1), extremesOfBlobs_ss(:, 1)), max(extremesOfBlobs(:, 2), extremesOfBlobs_ss(:, 2))];
                extremesOmIntTab_tot = [min(extremesOmIntTab(:, 1), extremesOmIntTab_ss(:, 1)), max(extremesOmIntTab(:, 2), extremesOmIntTab_ss(:, 2))];
            else
                lims_tot = lims;
                abs_lims_tot = abs_lims;

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

                extremesOfBlobs_tot = extremesOfBlobs;
                extremesOmIntTab_tot = extremesOmIntTab;
            end

            blobs = cell(1, tot_blobs);
            % Let's restrict blobs to the reached regions
            for ii = 1:tot_blobs
                blob = bl(ii).intm;
                blob = blob ./ sum(blob(:));

                if (overflow_b(ii) || overflow_ss_b(ii))
                    warning('Gt6DReconstructionAlgorithmFactory:BlobCreation', ...
                        'Blob: %03d --> Some orientations (%d, ss %d) will be missing', ...
                        ii, overflow_b(ii), overflow_ss_b(ii) )
                end

                blobs{ii} = blob(:, :, lims_tot(ii, 1):lims_tot(ii, 2));
                blobs{ii} = permute(blobs{ii}, [1 3 2]);
                blobs{ii} = cast(blobs{ii}, self.data_type);

                fprintf(' - ExremesBlobs %4d:%4d, ExtremesProjs %4d:%4d, Lims %3d:%3d, Size 1:%d, ', ...
                    extremesOfBlobs_tot(ii, :), extremesOmIntTab_tot(ii, :), lims_tot(ii, :), size(blobs{ii}, 3))

                intensExcl = sum(sum(sum(blob(:, :, [1:(lims_tot(ii, 1)-1), (lims_tot(ii, 2)+1):end]))));
                intensBlob = sum(sum(sum(blobs{ii})));
                fprintf('Blob: %03d -> Intensities: included %f, excluded %f\n', ...
                    ii, intensBlob, intensExcl);
            end

            [geometries, offsets] = Gt6DReconstructionAlgorithmFactory.makeBlobGeometries(grains_props, omega_int_tab, abs_lims_tot);

            if (~isempty(grains_ss))
                [geometries_ss, offsets_ss] = Gt6DReconstructionAlgorithmFactory.makeBlobGeometries(grains_props_ss, omega_int_tab_ss, abs_lims_tot);
            else
                geometries_ss = {};
                offsets_ss = {};
            end

            % Build Empty volumes
            blob_size = size(bl(1).intm);
            volume_size = Gt6DReconstructionAlgorithmFactory.getVolSize(blob_size, vol_type);

            num_geoms = numel(geometries);
            zero_volumes = arrayfun(@(x){zeros(volume_size, self.data_type)}, 1:num_geoms);

            algo = Gt6DBlobReconstructor(zero_volumes, blobs, ...
                'geometries', geometries, ...
                'offsets', offsets, ...
                'ss_geometries', geometries_ss, ...
                'ss_offsets', offsets_ss, ...
                'ss_coeffs', sampler.get_ss_coeffs() );
        end

        function [algo, blobs] = getSubBlobReconstructionAlgo(self, bl, selected_bls, sampler, vol_type, num_interp)
            fprintf('Extracting blobs on detector: ')
            c = tic();
            bl = bl(selected_bls);
            fprintf('\b\b (%f s).\n', toc(c));
            tot_blobs = numel(bl);

            [grains, grains_ss] = sampler.get_orientations();
            grains = [grains{:}];
            % Array of structures that contain all the info relative to
            % each orientation, for each blob
            grains_props = Gt6DReconstructionAlgorithmFactory.get_selected_props(grains, selected_bls);

            omega_int_tab = Gt6DReconstructionAlgorithmFactory.getOmegaIntegerTab(grains_props, sampler.parameters);

            [sub_blob_slices, proj_coeffs] = self.getSubBlobSlices(bl, num_interp);
            blobs = cell(size(sub_blob_slices));

            [slice_lims, lims, abs_lims, ~, extremesOfBlobs, extremesOmIntTab] = ...
                Gt6DReconstructionAlgorithmFactory.getSubBlobLims(omega_int_tab, bl, num_interp);

            [overflow_b, overflow_o] = Gt6DReconstructionAlgorithmFactory.compute_overflow(omega_int_tab, extremesOfBlobs);

            % We have to remove the orientations that don't project to any
            % blob
            class(grains_props)
            size(grains_props)
            good_or = overflow_o < tot_blobs;
            grains_props = grains_props(good_or);
            omega_int_tab = omega_int_tab(:, good_or);

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

                grains_props_ss = Gt6DReconstructionAlgorithmFactory.get_selected_props(grains_ss, selected_bls);
                omega_int_tab_ss = Gt6DReconstructionAlgorithmFactory.getOmegaIntegerTab(grains_props_ss, sampler.parameters);
                [slice_lims_ss, lims_ss, abs_lims_ss, ~, extremesOfBlobs_ss, extremesOmIntTab_ss] ...
                    = Gt6DReconstructionAlgorithmFactory.getSubBlobLims(omega_int_tab_ss, bl, num_interp);

                [overflow_ss_b, overflow_ss_o] = Gt6DReconstructionAlgorithmFactory.compute_overflow(omega_int_tab_ss, extremesOfBlobs_ss);

                good_ss_or = overflow_ss_o < tot_blobs;
                grains_props_ss = grains_props_ss(good_ss_or);
                omega_int_tab_ss = omega_int_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))];

                extremesOfBlobs_tot = [min(extremesOfBlobs(:, 1), extremesOfBlobs_ss(:, 1)), max(extremesOfBlobs(:, 2), extremesOfBlobs_ss(:, 2))];
                extremesOmIntTab_tot = [min(extremesOmIntTab(:, 1), extremesOmIntTab_ss(:, 1)), max(extremesOmIntTab(:, 2), extremesOmIntTab_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));

                extremesOfBlobs_tot = extremesOfBlobs;
                extremesOmIntTab_tot = extremesOmIntTab;
            end

            % Let's restrict blobs to the reached regions
            for ii = 1:tot_blobs
                if (overflow_b(ii) || overflow_ss_b(ii))
                    warning('Gt6DReconstructionAlgorithmFactory:BlobCreation', ...
                        'Blob: %03d --> Some orientations (%d, ss %d) will be missing', ...
                        ii, overflow_b(ii), overflow_ss_b(ii) )
                end

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

                fprintf(' - ExremesBlobs %4d:%4d, ExtremesProjs %4d:%4d, Lims %3d:%3d, SliceLims %d:%d, Size 1:%d, ', ...
                    extremesOfBlobs_tot(ii, :), extremesOmIntTab_tot(ii, :), lims_tot(ii, :), slice_lims_tot(ii, :), size(sub_blob_slices{ii}, 3))

                blobs{ii} = sub_blob_slices{ii}(:, :, slice_lims_tot(ii, 1):slice_lims_tot(ii, 2));
                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('Blob: %03d -> Intensities: included %f, excluded %f\n', ...
                    ii, intensIncl, intensExcl);
            end

            [geometries, offsets] = Gt6DReconstructionAlgorithmFactory.makeSubBlobGeometries( ...
                grains_props, proj_coeffs, extremesOfBlobs_tot, omega_int_tab, abs_lims_tot, slice_lims_tot);

            if (~isempty(grains_ss))
                [geometries_ss, offsets_ss] = Gt6DReconstructionAlgorithmFactory.makeSubBlobGeometries( ...
                    grains_props_ss, proj_coeffs, extremesOfBlobs_tot, omega_int_tab_ss, abs_lims_tot, slice_lims_tot);
            else
                geometries_ss = {};
                offsets_ss = {};
            end

            % Build Empty volumes
            blob_size = size(bl(1).intm);
            volume_size = Gt6DReconstructionAlgorithmFactory.getVolSize(blob_size, vol_type);

            num_geoms = numel(geometries);
            zero_volumes = arrayfun(@(x){zeros(volume_size, self.data_type)}, 1:num_geoms);

            algo = Gt6DBlobReconstructor(zero_volumes, blobs, ...
                'geometries', geometries, ...
                'offsets', offsets, ...
                'ss_geometries', geometries_ss, ...
                'ss_offsets', offsets_ss, ...
                'ss_coeffs', sampler.get_ss_coeffs() );
        end
    end

    methods (Access = protected)
        function [subBlobSlices, projCoeffs] = getSubBlobSlices(self, blobs, numInterp)
            subBlobSlices = cell(size(blobs));
            subBlobCoeffs = cell(size(blobs));
            projCoeffs = cell(size(blobs));

            blocky_groups = false;

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

                [blobSize(1), blobSize(2), blobSize(3)] = size(blobs(ii).intm);
                % Number of interpolated slices should be at least enough
                % to accomodate all the slices of the blob
                numSlices = ceil((blobSize(3) -1) / numInterp) +1;

                covered_slices_size = (numSlices-1) * numInterp + 1;

%                 num_chars = num_chars + fprintf(' Blob size %d, numSlices %d, convertedSliceSize %d', ...
%                     blobSize(3), numSlices, covered_slices_size);

                subBlobSlices{ii} = zeros(blobSize(1), blobSize(2), numSlices, self.data_type);
                for n = numSlices:-1:1
                    coeffs = abs( (numInterp * (n-1) +1) - (1:covered_slices_size));
                    coeffs = 1 - coeffs / numInterp;

                    % these will be the blob slices that relate to the
                    % given ii interpolated slice
                    if (blocky_groups)
                        goodCoeffs = find(coeffs >= 0.5);
                    else
                        goodCoeffs = find(coeffs > 0);
                    end
                    % The corresponding coefficients
                    coeffs = coeffs(goodCoeffs);

                    if (blocky_groups)
                        coeffs(1:end) = 1;
                        if (mod(numInterp, 2) == 0)
                            coeffs(1) = 0.5;
                            coeffs(end) = 0.5;
                        end
                    end

                    safe_coeffs_pos = goodCoeffs <= blobSize(3);
                    safe_good_coeffs = goodCoeffs(safe_coeffs_pos);
                    safe_coeffs = coeffs(safe_coeffs_pos);
                    % Let's build the interpolated slices
                    blob = blobs(ii).intm;
                    blob = blob ./ sum(blob(:));
                    for w = 1:numel(safe_coeffs)
                        subBlobSlices{ii}(:, :, n) = ...
                            subBlobSlices{ii}(:, :, n) + safe_coeffs(w) .* blob(:, :, safe_good_coeffs(w));
                    end
                    % Save the metadata
                    subBlobCoeffs{ii}(n) = struct('indexes', goodCoeffs, 'coeffs', coeffs);
                end

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

                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 subBlobCoeffs is the other way
            % around.
            for ii = 1:numel(subBlobSlices)
                numSlices = numel(subBlobCoeffs{ii});

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

    methods (Access = protected, Static)
        function grains_props = get_selected_props(grains, blobs_on_det)
%             grains_props = [grains(:).allblobs];

            for ii = numel(grains):-1:1
                all_props = grains(ii).allblobs;
                grains_props(ii) = Gt6DReconstructionAlgorithmFactory.selectBlobsFromGrain(...
                    all_props, blobs_on_det);
            end
        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 volume_size = getVolSize(blob_size, vol_type)
%             volume_size = floor(sqrt(blob_size(1:2) * blob_size(1:2)') /2)+1;
            volume_size = [blob_size(1), blob_size(1:2)];
%             volume_size = volume_size * 0.7;
            volume_size = volume_size * 0.8;
%             volume_size = volume_size * 1.2;
            volume_size = volume_size - mod(volume_size, 8);

            switch (vol_type)
                case '3D'
%                     volume_size = [volume_size, volume_size, volume_size];
                case '2D'
%                     volume_size = [volume_size, volume_size, 8];
                    volume_size(3) = 8;
                otherwise
                    error('Gt6DReconstructionAlgorithmFactory:getVolSize:wrong_argument', ...
                        'Volume type should be either "2D" or "3D"')
            end
        end

        function omegaIntegerTab = getOmegaIntegerTab(grainsProps, parameters)
            omega_step = 180 / parameters.acq.nproj;
            % Tabular of omegas: horizontal direction the n geometries,
            % vertical direction the omegas relative to each blob
            omegaTab = [grainsProps(:).omega];
            omegaSlicesTab = omegaTab / omega_step;
            % Images numbers
            omegaIntegerTab = round(omegaSlicesTab);
        end

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

        function [lims, absLims, overflow, extremesOfBlobs, extremesOmIntTab] = getBlobLims(omegaIntegerTab, bl)
            % Let's add extreme values
            extremesOmIntTab = [ ...
                min(omegaIntegerTab, [], 2), ...
                max(omegaIntegerTab, [], 2) ];

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

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

            overflow = ...
                (extremesOmIntTab(:, 1) < extremesOfBlobs(:, 1)) ...
                | (extremesOmIntTab(:, 2) > extremesOfBlobs(:, 2));
        end

        function [sliceLims, lims, absLims, overflow, extremesOfBlobs, extremesOmIntTab] = getSubBlobLims(omegaIntegerTab, bl, numInterp)
            [~, ~, ~, extremesOfBlobs, extremesOmIntTab] = Gt6DReconstructionAlgorithmFactory.getBlobLims(omegaIntegerTab, bl);

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

            extremesOfBlobs(:, 2) = extremesOfBlobs(:, 1) + (blobs_sizes-1) ...
                + mod(numInterp - mod((blobs_sizes-1), numInterp), numInterp);

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

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

            overflow = ...
                (extremesOmIntTab(:, 1) < extremesOfBlobs(:, 1)) ...
                | (extremesOmIntTab(:, 2) > extremesOfBlobs(:, 2));
        end

        function [omega_int_tab, omega_int_coeffs] = getOmegaIntegersAndCoeffs(grainsProps, parameters)
            omega_step = 180 / parameters.acq.nproj;
            % Tabular of omegas: horizontal direction the n geometries,
            % vertical direction the omegas relative to each blob
            omega_tab = [grainsProps(:).omega];
            omega_slices_tab = omega_tab / omega_step;
            % Images numbers
            omega_int_tab(:, :, 1) = floor(omega_slices_tab);
            omega_int_coeffs(:, :, 1) = 1 - (omega_slices_tab - omega_integer_min_tab);
            omega_int_tab(:, :, 2) = ceil(omega_slices_tab);
            omega_int_coeffs(:, :, 2) = 1 - omega_int_coeffs(:, :, 1);
        end

        function [geometries, offsets] = makeBlobGeometries(grains_props, omega_int_tab, abs_lims)
            % ASTRA Geometry
            geometries = {grains_props(:).proj_geom};
            num_geoms = numel(geometries);
            offsets = cell(1, num_geoms);

            % Compute Offsets
            rescaled_omegas = omega_int_tab - repmat(abs_lims(:, 1), [1 num_geoms]) + 1;
            exceding_omegas = omega_int_tab - repmat(abs_lims(:, 2), [1 num_geoms]);

            % Assign offsets
            for ii = 1:num_geoms
                valid_proj_dirs = ...
                    (rescaled_omegas(:, ii) >= 1) ...
                    & (exceding_omegas(:, ii) < 1);

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

                offsets{ii} = struct('blob_offsets', find(valid_proj_dirs), ...
                    'proj_offsets', rescaled_omegas(valid_proj_dirs, ii), ...
                    'proj_coeffs', ones(numel(find(valid_proj_dirs)), 1), ...
                    'sino_offsets', 1:numel(find(valid_proj_dirs)) );

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

        function [geometries, offsets] = makeSubBlobGeometries( ...
                grains_props, proj_coeffs, extremesOfBlobs, omega_int_tab, abs_lims, slice_lims)
            % ASTRA Geometry
            geometries = {grains_props(:).proj_geom};
            num_geoms = numel(geometries);

%             extremesOfBlobs = vertcat( bl(:).bbwim );
            slice_pos_in_blob = omega_int_tab - repmat(extremesOfBlobs(:, 1), [1 num_geoms]) + 1;

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

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

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

                valid_proj_indxes = find(valid_proj_dirs)';
                for jj = 1:numel(valid_proj_indxes)
                    n = valid_proj_indxes(jj);

                    slices = [proj_coeffs{n}(slice_pos_in_blob(n, ii)).slices];
                    coeffs = [proj_coeffs{n}(slice_pos_in_blob(n, ii)).coeffs];

                    slices = slices - slice_lims(n, 1) +1;

                    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
    end
end