PSF: small fix for non-symmetric PSFs

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

5_reconstruction/gtAstra3DTV.m

@@ -79,7 +79,7 @@ function [volume, res_norm] = gtAstra3DTV(proj, varargin)
         psf_obj.set_psf_direct(proj.psf, [proj.num_cols, proj.num_rows]);
 
         fp = @(x)psf_obj.apply_psf_direct(astra_mex_direct_c('FP3D', astra_projector_id, x));
-        bp = @(x)astra_mex_direct_c('BP3D', astra_projector_id, psf_obj.apply_psf_transpose(x));
+        bp = @(x)astra_mex_direct_c('BP3D', astra_projector_id, psf_obj.apply_psf_adjoint(x));
     else
         fp = @(x)astra_mex_direct_c('FP3D', astra_projector_id, x);
         bp = @(x)astra_mex_direct_c('BP3D', astra_projector_id, x);

zUtil_Imaging/GtPSF.m

 classdef GtPSF < handle
     properties
         psf_direct = [];
-        psf_transpose = [];
+        psf_adjoint = [];
 
         otf_direct = [];
-        otf_transpose = [];
+        otf_adjoint = [];
 
         is_symmetric = true;
         is_projection_invariant = true;
@@ -24,38 +24,18 @@ classdef GtPSF < handle
         function set_psf_direct(self, psf_d, image_size)
             self.reset();
 
-            % let's be sure that it will be in ASTRA's projection
-            % convention
-            if (size(psf_d, 1) == size(psf_d, 2))
-                psf_d = permute(psf_d, [1 3 2]);
-            elseif (size(psf_d, 1) ~= size(psf_d, 3))
-                error([mfilename ':wrong_argument'], ...
-                    'PSFs should be in the form of _square images_, with odd edges')
-            elseif (mod(size(psf_d, 1), 2) ~= 1)
-                error([mfilename ':wrong_argument'], ...
-                    'PSFs should be in the form of square images, with _odd edges_')
-            end
+            [psf_d, psf_a, self.is_symmetric, self.is_projection_invariant] ...
+                = self.check_single_psf(psf_d);
 
-            % let's renormalize the PSFs
-            psfs_norm = 1 ./ sum(sum(abs(psf_d), 1), 3);
-            psf_d = bsxfun(@times, psf_d, psfs_norm);
-
-            self.psf_direct = cast(psf_d, self.data_type);
-            if (exist('image_size', 'var') && ~isempty(image_size))
-                self.image_size = image_size;
-            end
-
-            % Let's find out whether they are symmetric or not
-            psf_t = permute(psf_d, [3 2 1]);
-            self.is_symmetric = all(abs(psf_d(:) - psf_t(:)) < eps('single'));
+            self.psf_direct = psf_d;
             if (~self.is_symmetric)
-                self.psf_transpose = psf_t;
+                self.psf_adjoint = psf_a;
             end
 
-            self.is_projection_invariant = size(psf_d, 2) == 1;
+            if (exist('image_size', 'var') && ~isempty(image_size))
+                self.image_size = image_size;
 
-            % if we know the images sizes, we can already compute the OTF
-            if (~isempty(self.image_size))
+                % if we know the images sizes, we can already compute the OTF
                 self.init_otfs();
             end
         end
@@ -70,7 +50,7 @@ classdef GtPSF < handle
             end
         end
 
-        function imgs = apply_psf_transpose(self, imgs)
+        function imgs = apply_psf_adjoint(self, imgs)
             self.check_incoming_images(imgs);
 
             if (self.use_otf)
@@ -81,17 +61,44 @@ classdef GtPSF < handle
         end
     end
 
-    methods (Access = private)
+    methods (Access = protected)
         function reset(self)
             self.psf_direct = [];
-            self.psf_transpose = [];
+            self.psf_adjoint = [];
 
             self.otf_direct = [];
-            self.otf_transpose = [];
+            self.otf_adjoint = [];
 
             self.image_size = [];
         end
 
+        function [psf_d, psf_a, is_sym, is_p_inv] = check_single_psf(self, psf_d)
+            % let's be sure that it will be in ASTRA's projection
+            % convention
+            if (size(psf_d, 1) == size(psf_d, 2))
+                psf_d = permute(psf_d, [1 3 2]);
+            elseif (size(psf_d, 1) ~= size(psf_d, 3))
+                error([mfilename ':wrong_argument'], ...
+                    'PSFs should be in the form of _square images_, with odd edges')
+            elseif (mod(size(psf_d, 1), 2) ~= 1)
+                error([mfilename ':wrong_argument'], ...
+                    'PSFs should be in the form of square images, with _odd edges_')
+            end
+
+            % let's renormalize the PSFs
+            psfs_norm = 1 ./ sum(sum(abs(psf_d), 1), 3);
+            psf_d = bsxfun(@times, psf_d, psfs_norm);
+
+            psf_d = cast(psf_d, self.data_type);
+
+            % Let's find out whether they are symmetric or not
+            psf_a = flip(flip(psf_d, 1), 3);
+
+            is_sym = all(abs(psf_d(:) - psf_a(:)) < eps('single'));
+
+            is_p_inv = size(psf_d, 2) == 1;
+        end
+
         function otf = init_single_otf(self, psf)
             num_psfs = size(psf, 2);
             psf_edge = (size(psf, 1) - 1) / 2;
@@ -122,7 +129,7 @@ classdef GtPSF < handle
             self.otf_direct = self.init_single_otf(self.psf_direct);
 
             if (~self.is_symmetric)
-                self.otf_transpose = self.init_single_otf(self.psf_transpose);
+                self.otf_adjoint = self.init_single_otf(self.psf_adjoint);
             end
         end
 
@@ -137,7 +144,7 @@ classdef GtPSF < handle
             imgs_size = [size(imgs, 1), size(imgs, 3)];
             if (isempty(self.image_size))
                 self.otf_direct = [];
-                self.otf_transpose = [];
+                self.otf_adjoint = [];
 
                 self.image_size = imgs_size;
             elseif (any(self.image_size ~= imgs_size))
@@ -146,7 +153,7 @@ classdef GtPSF < handle
                     self.image_size, imgs_size)
 
                 self.otf_direct = [];
-                self.otf_transpose = [];
+                self.otf_adjoint = [];
 
                 self.image_size = imgs_size;
             end
@@ -171,7 +178,7 @@ classdef GtPSF < handle
             if (is_direct || self.is_symmetric)
                 imgs = bsxfun(@times, imgs, self.otf_direct);
             else
-                imgs = bsxfun(@times, imgs, self.otf_transpose);
+                imgs = bsxfun(@times, imgs, self.otf_adjoint);
             end
 
             imgs = ifftshift(imgs, 1);
@@ -193,7 +200,7 @@ classdef GtPSF < handle
             if (is_direct || self.is_symmetric)
                 psf = self.psf_direct;
             else
-                psf = self.psf_transpose;
+                psf = self.psf_adjoint;
             end
 
             if (self.is_projection_invariant)