6D-reconstruction: better handling of PSFs in reconstructions

This allows to define a global PSF for reconstruction, but it also handles the case where different grains (in cluster reconstructions) were to have different PSFs for some reason.
Ideally every blob should have a different PSF.

Signed-off-by: Nicola Vigano <nicola.vigano@esrf.fr>

zUtil_Deformation/Gt6DReconstructionAlgorithmFactory.m

@@ -57,6 +57,10 @@ classdef Gt6DReconstructionAlgorithmFactory < handle
             blobs = algo_params.blobs;
             or_groups = [1, sampler.get_total_num_orientations()];
 
+            if (isempty(algo_params.psf) && isfield(self.parameters.rec.grains.options, 'psf'))
+                algo_params.psf = { permute(self.parameters.rec.grains.options.psf, [1 3 2]) };
+            end
+
             algo = Gt6DBlobReconstructor(volume_size, blobs, ...
                 'data_type', self.data_type, ...
                 'geometries', algo_params.geometries, ...
@@ -128,10 +132,47 @@ classdef Gt6DReconstructionAlgorithmFactory < handle
 
             blobs = cat(1, algo_params(:).blobs);
 
-            psf = cat(1, algo_params(:).psf);
-%             psf = fspecial('gaussian', 21, 0.75) + 2 * fspecial('gaussian', 21, 4);
-%             psf = psf / sum(psf(:));
-%             psf = { permute(psf, [1 3 2]) };
+            % Handling PSFs
+            num_psfs = arrayfun(@(x)numel(x.psf), algo_params);
+            num_bls = arrayfun(@(x)numel(x.blobs), algo_params);
+            if (all(num_psfs == 0))
+                % If none of the regions has a dedicated PSF-stack, we
+                % check in the options
+                if (isfield(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 all orientation-space regions\n\n')
+                    psf = {};
+                end
+            else
+                psf = {};
+                fprintf('\nDifferent PSFs used for the different orientation-space regions:\n')
+                for ii_o_reg = 1:num_regions
+                    if (num_psfs(ii_o_reg) > 1)
+                        fprintf('%2d) PSFs: %d, Blobs: %d\n', ii_o_reg, num_psfs(ii_o_reg), num_bls(ii_o_reg))
+                        if (num_psfs(ii_o_reg) ~= num_bls(ii_o_reg))
+                            error('Gt6DReconstructionAlgorithmFactory:wrong_psf_settings', ...
+                                'Orientation-space region %d has a different number of PSFs from the number of blobs', ...
+                                ii_o_reg);
+                        end
+                        new_psfs = algo_params(ii_o_reg).psf;
+                    else
+                        fprintf('%2d) PSFs: %d -> %d, Blobs: %d\n', ii_o_reg, num_psfs(ii_o_reg), num_bls(ii_o_reg), num_bls(ii_o_reg))
+                        if (num_psfs(ii_o_reg) == 0)
+                            if (isfield(self.parameters.rec.grains.options, 'psf'))
+                                algo_params(ii_o_reg).psf = { permute(self.parameters.rec.grains.options.psf, [1 3 2]) };
+                            else
+                                algo_params(ii_o_reg).psf = { ones(1, 1, 1) };
+                            end
+                        end
+                        ones_exp = ones(num_bls(ii_o_reg), 1);
+                        new_psfs = algo_params(ii_o_reg).psf(ones_exp);
+                    end
+                    psf = [psf, new_psfs]; %#ok<AGROW>
+                end
+                fprintf('\n')
+            end
 
             geometries = cat(1, algo_params(:).geometries);
             geometries_ss = cat(1, algo_params(:).geometries_ss);