GtGrainODFuvwPatchSolver.m

classdef GtGrainODFuvwPatchSolver < GtGrainODFuvwSolver
    properties
        shape_functions_extracted = {};
    end

    methods (Access = public)
        function self = GtGrainODFuvwPatchSolver(parameters, varargin)
            self = self@GtGrainODFuvwSolver(parameters, varargin{:});
        end
    end

    methods (Access = public)
        function build_sinogram(self)
            self.build_sinogram@GtGrainODFuvwSolver()

            self.sino = reshape(self.sino, self.size_sino);
        end

        function build_orientation_sampling(self, ref_gr, det_index, oversize, oversampling, shape_functions)
            self.build_orientation_sampling@GtGrainODFuvwSolver(ref_gr, det_index, oversize, oversampling, shape_functions);

            for ii_o = 1:self.sampler.get_total_num_orientations()
                for ii_v = 1:size(self.voxel_centers, 1)
                    self.shape_functions_extracted{ii_o, ii_v} = {self.shape_functions{ii_v}{ii_o}(:).intm};
                end
            end
        end

        function build_projection_matrices_sf_uvw(self)
            ref_sel = self.sampler.selected;

            grid_gr = self.sampler.get_orientations();
            num_orients = numel(grid_gr);

            num_voxels = size(self.voxel_centers, 1);
            self.S = cell(num_voxels, 1);
            self.St = cell(num_voxels, 1);

            bls = self.sampler.ref_gr.proj(self.sampler.detector_index).bl(ref_sel);

            bls_uvw_orig = cat(1, bls(:).bbuim, bls(:).bbvim, bls(:).bbwim);
            bls_uvw_orig = reshape(bls_uvw_orig, [], 3, 2);
            bls_uvw_orig = bls_uvw_orig(:, :, 1);
            bls_uvw_orig(:, 3) = bls_uvw_orig(:, 3) - self.pre_paddings + 1;

            for ii_v = 1:num_voxels
                num_chars = fprintf('%03d/%03d', ii_v, num_voxels);

                sf_uvw_shift = cell(num_orients, 1);
                for ii_o = 1:num_orients
                    sf_uvw_shift{ii_o} = cat(1, ...
                        self.shape_functions{ii_v}{ii_o}(:).bbuim, ...
                        self.shape_functions{ii_v}{ii_o}(:).bbvim, ...
                        self.shape_functions{ii_v}{ii_o}(:).bbwim ...
                        );
                    sf_uvw_shift{ii_o} = reshape(sf_uvw_shift{ii_o}, [], 3, 2);
                    sf_uvw_shift{ii_o} = sf_uvw_shift{ii_o}(:, :, 1) - bls_uvw_orig;
                end
                sf_uvw_shift = cat(3, sf_uvw_shift{:});

                self.S{ii_v} = sf_uvw_shift;

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

%         function build_projection_matrices_sf_uvw_approx(self)
%             ref_sel = self.sampler.selected;
% 
%             det_ind = self.sampler.detector_index;
% 
%             detgeo = self.parameters.detgeo(det_ind);
%             labgeo = self.parameters.labgeo;
%             samgeo = self.parameters.samgeo;
% 
%             grid_gr = self.sampler.get_orientations();
%             num_orients = numel(grid_gr);
% 
%             num_voxels = size(self.voxel_centers, 1);
%             self.S = cell(num_voxels, 1);
%             self.St = cell(num_voxels, 1);
% 
%             % Preparing for computing the UV deviations due to the XYZ
%             % positions of the real-space voxels
%             rot_l2s = cell(num_orients, 1);
%             for ii_o = 1:num_orients
%                 rot_l2s{ii_o} = grid_gr{ii_o}.allblobs(det_ind).srot(:, :, ref_sel);
%             end
%             rot_l2s = cat(3, rot_l2s{:});
%             rot_s2l = permute(rot_l2s, [2 1 3]);
% 
%             dveclab = cell(num_orients, 1);
%             for ii_o = 1:num_orients
%                 dveclab{ii_o} = grid_gr{ii_o}.allblobs(det_ind).dvec(ref_sel, :);
%             end
%             dveclab = cat(1, dveclab{:});
% 
%             uv_orig = cell(num_orients, 1);
%             for ii_o = 1:num_orients
%                 uv_orig{ii_o} = grid_gr{ii_o}.allblobs(det_ind).detector.uvw(ref_sel, 1:2);
%             end
%             uv_orig = cat(3, uv_orig{:});
%
%             bls = self.sampler.ref_gr.proj(self.sampler.detector_index).bl(ref_sel);
%
%             bls_uvw_orig = cat(1, bls(:).bbuim, bls(:).bbvim, bls(:).bbwim);
%             bls_uvw_orig = reshape(bls_uvw_orig, [], 3, 2);
%             bls_uvw_orig = bls_uvw_orig(:, :, 1);
%             bls_uvw_orig(:, 3) = bls_uvw_orig(:, 3) - self.pre_paddings + 1;
% 
%             sf_uvw_shift = cell(num_orients, 1);
%             for ii_o = 1:num_orients
%                 sf_uvw_shift{ii_o} = cat(1, ...
%                     self.shape_functions{ii_o}(:).bbuim, ...
%                     self.shape_functions{ii_o}(:).bbvim, ...
%                     self.shape_functions{ii_o}(:).bbwim ...
%                     );
%                 sf_uvw_shift{ii_o} = reshape(sf_uvw_shift{ii_o}, [], 3, 2);
%                 sf_uvw_shift{ii_o} = sf_uvw_shift{ii_o}(:, :, 1) - bls_uvw_orig;
%             end
%             sf_uvw_shift = cat(3, sf_uvw_shift{:});
% 
%             for ii_v = 1:num_voxels
%                 num_chars = fprintf('%03d/%03d', ii_v, num_voxels);
% 
%                 gcsam_v = self.voxel_centers(ii_v, :);
%                 gcsam_v = gcsam_v(ones(1, size(rot_s2l, 3)), :);
%                 gclab_v = gtGeoSam2Lab(gcsam_v, rot_s2l, labgeo, samgeo, false);
% 
%                 uv_tot = gtFedPredictUVMultiple([], dveclab', gclab_v', ...
%                     detgeo.detrefpos', detgeo.detnorm', detgeo.Qdet, ...
%                     [detgeo.detrefu, detgeo.detrefv]');
%                 uv_tot = reshape(uv_tot, 2, [], num_orients);
% 
%                 uv_shift = permute(uv_tot, [2 1 3]) - uv_orig;
%                 self.S{ii_v} = sf_uvw_shift + [uv_shift, zeros(size(uv_shift, 1), 1, num_orients)];
% 
%                 fprintf(repmat('\b', [1 num_chars]));
%             end
%         end
    end

    methods (Access = protected)
        function y = fp(self, x)
            num_blobs = self.get_num_blobs();
            num_orients = size(self.shape_functions_extracted, 1);
            num_voxels = size(self.shape_functions_extracted, 2);

            y = gtMathsGetSameSizeZeros(self.sino, self.data_type);

            for ii_v = 1:num_voxels
                uvw_shifts = self.S{ii_v};

                for ii_o = 1:num_orients
                    bls = self.shape_functions_extracted{ii_o, ii_v};

                    for ii_b = 1:num_blobs
                        bls{ii_b} = bls{ii_b} .* x(ii_o, ii_v);
                    end
                    y = gtPlaceSubVolumes(y, bls, [uvw_shifts(:, :, ii_o), (0:num_blobs-1)']);
                end
            end
        end

        function y = bp(self, x)
            num_blobs = self.get_num_blobs();
            num_orients = size(self.shape_functions_extracted, 1);
            num_voxels = size(self.shape_functions_extracted, 2);

            y = zeros([num_orients, num_voxels], self.data_type);

            for ii_v = 1:num_voxels
                uvw_shifts = self.S{ii_v};

                for ii_o = 1:num_orients
                    bls = self.shape_functions_extracted{ii_o, ii_v};

                    mask_sino = zeros(self.size_sino, self.data_type);
                    mask_sino = gtPlaceSubVolumes(mask_sino, bls, [uvw_shifts(:, :, ii_o), (0:num_blobs-1)']);

                    mask_sino = mask_sino .* x;

                    y(ii_o, ii_v) = gtMathsSumNDVol(mask_sino);
                end
            end
        end

%         function y = fp_approx(self, x)
%             num_blobs = self.get_num_blobs();
% 
%             y = gtMathsGetSameSizeZeros(self.sino);
% 
%             for ii_o = 1:numel(self.shape_functions)
%                 sfs = self.shape_functions{ii_o};
%                 for ii_b = 1:num_blobs
%                     sf = sfs(ii_b).intm;
%                     for ii_v = 1:numel(self.S)
%                         sf_rescaled = sf .* x(ii_o, ii_v);
%                         y = gtPlaceSubVolume(y, sf_rescaled, [self.S{ii_v}(ii_b, :, ii_o), ii_b], [], 'summed');
%                     end
%                 end
%             end
%         end
% 
%         function y = bp_approx(self, x)
%             num_blobs = self.get_num_blobs();
% 
%             num_orients = numel(self.shape_functions);
%             num_voxels = numel(self.S);
% 
%             y = zeros([num_orients, numel(self.S)]);
% 
%             for ii_v = 1:num_voxels
%                 for ii_o = 1:num_orients
%                     sfs = self.shape_functions{ii_o};
% 
%                     mask_sino = gtMathsGetSameSizeZeros(self.sino);
% 
%                     for ii_b = 1:num_blobs
%                         sf = sfs(ii_b).intm;
%                         mask_sino = gtPlaceSubVolume(mask_sino, sf, [self.S{ii_v}(ii_b, :, ii_o), ii_b], [], 'summed');
%                     end
% 
%                     mask_sino = mask_sino .* x;
%                     y(ii_o, ii_v) = sum(mask_sino(:));
%                 end
%             end
%         end
    end
end