PSF: little improvements to the PSF producing functions

Signed-off-by: Nicola VIGANO <N.R.Vigano@cwi.nl>

zUtil_Imaging/gtESF2LSF.m

@@ -6,33 +6,41 @@ function lsf_line = gtESF2LSF(img, varargin)
         'keep_oversampling', true, 'verbose', true);
     conf = parse_pv_pairs(conf, varargin);
 
-    p = gtLoadParameters();
-
     if (ischar(img))
         img = GtVolView.loadVolume(img);
     elseif (numel(img) == 1)
+        p = gtLoadParameters();
+
         filename = sprintf('refHST%04d.edf', img);
         filename = fullfile(p.acq.dir, '0_rawdata', p.acq.name, filename);
         img = edf_read(filename);
     end
 
     if (isempty(conf.bb))
-        bb_half_size = ceil(p.acq.bbdir(3:4) / 2);
-        bb = [p.acq.bbdir(1:2), bb_half_size];
-        if (isempty(conf.direction) || strcmpi(conf.direction, 'top-bottom'))
-            bb(1) = bb(1) + ceil(bb_half_size(1) / 2);
-            bb(2) = bb(2) - ceil(bb_half_size(2) / 2);
-        elseif (strcmpi(conf.direction, 'left-right'))
-            bb(1) = bb(1) - ceil(bb_half_size(1) / 2);
-            bb(2) = bb(2) + ceil(bb_half_size(2) / 2);
-        elseif (strcmpi(conf.direction, 'right-left'))
-            bb(1) = bb(1) + 3 * ceil(bb_half_size(1) / 2);
-            bb(2) = bb(2) + ceil(bb_half_size(2) / 2);
-        else
-            bb(1) = bb(1) + ceil(bb_half_size(1) / 2);
-            bb(2) = bb(2) + 3 * ceil(bb_half_size(2) / 2);
+        try
+            if (~exist('p', 'var') || isempty(p))
+                p = gtLoadParameters();
+            end
+
+            bb_half_size = ceil(p.acq.bbdir(3:4) / 2);
+            bb = [p.acq.bbdir(1:2), bb_half_size];
+            if (isempty(conf.direction) || strcmpi(conf.direction, 'top-bottom'))
+                bb(1) = bb(1) + ceil(bb_half_size(1) / 2);
+                bb(2) = bb(2) - ceil(bb_half_size(2) / 2);
+            elseif (strcmpi(conf.direction, 'left-right'))
+                bb(1) = bb(1) - ceil(bb_half_size(1) / 2);
+                bb(2) = bb(2) + ceil(bb_half_size(2) / 2);
+            elseif (strcmpi(conf.direction, 'right-left'))
+                bb(1) = bb(1) + 3 * ceil(bb_half_size(1) / 2);
+                bb(2) = bb(2) + ceil(bb_half_size(2) / 2);
+            else
+                bb(1) = bb(1) + ceil(bb_half_size(1) / 2);
+                bb(2) = bb(2) + 3 * ceil(bb_half_size(2) / 2);
+            end
+            bb(3:4) = bb(3:4) + 1 - mod(bb(3:4), 2);
+        catch
+            bb = [0, 0, size(img, 2), size(img, 1)];
         end
-        bb(3:4) = bb(3:4) + 1 - mod(bb(3:4), 2);
 
         conf.bb = get_roi_bb(img, 'ESF Bounding Box', bb);
     end

zUtil_Imaging/gtLSFs2PSF.m

 function psf = gtLSFs2PSF(lsfs, varargin)
-    conf = struct('method', 'fourier');
+    conf = struct('method', 'fourier', 'downsampling', 1);
     conf = parse_pv_pairs(conf, varargin);
 
     num_lsfs = numel(lsfs);
 
+    data = cat(1, lsfs.data);
+    num_data_pixels = size(data, 2);
+
     if (num_lsfs == 1)
         % We're assuming here that the only LSF given is the radial
         % component of a radially symmetric PSF
-        data = reshape(lsfs.data, [], 1);
+        data = reshape(data, [], 1);
+        data = (data + flipud(data)) / 2;
+
         switch (lower(conf.method))
             case 'astra'
             case 'iradon'
                 data = data(:, ones(180, 1));
                 psf_inner = iradon(data, 1:180,'linear', 'Shepp-Logan');
-                psf_upsize = [numel(lsfs.data), numel(lsfs.data)];
+                psf_upsize = [num_data_pixels, num_data_pixels];
                 psf = zeros(psf_upsize, lsfs.data_type);
                 shifts = (psf_upsize - 1) / 2 - size(psf_inner) / 2 + 1;
                 psf(1+shifts(1):size(psf_inner, 1)+shifts(1), 1+shifts(2):size(psf_inner, 2)+shifts(2)) = psf_inner;
@@ -22,7 +27,7 @@ function psf = gtLSFs2PSF(lsfs, varargin)
                 data_f = fft(data_f);
                 data_f = fftshift(data_f);
 
-                psf_edge = (numel(lsfs.data) - 1) / 2;
+                psf_edge = (num_data_pixels - 1) / 2;
                 [xx, yy] = ndgrid(-psf_edge:psf_edge, -psf_edge:psf_edge);
                 rr = sqrt(xx .^ 2 + yy .^ 2);
                 psf_f = interp1(-psf_edge:psf_edge, data_f, rr, 'spline', 0);
@@ -50,13 +55,11 @@ function psf = gtLSFs2PSF(lsfs, varargin)
         end
 
         % The LSFs given here, are Fourier samplings of the 2D PSF
-        data = cat(1, lsfs.data);
-
         data_f = ifftshift(data, 2);
         data_f = fft(data_f, [], 2);
         data_f = fftshift(data_f, 2);
 
-        psf_edge = (size(data, 2) - 1) / 2;
+        psf_edge = (num_data_pixels - 1) / 2;
 
         dirs = cat(1, lsfs.direction);
         angles = mod(atan2d(dirs(:, 2), dirs(:, 1)), 360);
@@ -87,9 +90,30 @@ function psf = gtLSFs2PSF(lsfs, varargin)
         psf = real(fftshift(psf));
     end
 
-    if (lsfs(1).keep_oversampling)
-        psf = reshape(psf, lsfs(1).oversampling, lsfs(1).size, lsfs(1).oversampling, lsfs(1).size);
-        psf = squeeze(sum(sum(psf, 1), 3));
+    if (lsfs(1).keep_oversampling) % It was kept before, so now we get rid of it
+        to_be_summed = conf.downsampling * lsfs(1).oversampling;
+    else
+        to_be_summed = conf.downsampling;
+    end
+    extra_pixels = mod(to_be_summed - mod(num_data_pixels, to_be_summed), to_be_summed);
+    num_data_pixels = num_data_pixels + extra_pixels;
+
+    if (mod(conf.downsampling, 2))
+        if (mod(lsfs(1).size, conf.downsampling))
+            padsize = (extra_pixels - 1) / 2;
+            psf = padarray(psf, [padsize, padsize], 'both');
+        end
+    else
+        half_in_size = (size(psf) - 1) / 2;
+        half_out_size = (num_data_pixels - 1) / 2;
+        [xx, yy] = ndgrid(-half_out_size:half_out_size, -half_out_size:half_out_size);
+        psf = interp2(-half_in_size(1):half_in_size(1), -half_in_size(2):half_in_size(2), psf, xx, yy, 'spline', 0);
+        size(psf)
     end
+
+    final_psf_size = ceil(lsfs(1).size / conf.downsampling);
+
+    psf = reshape(psf, to_be_summed, final_psf_size, to_be_summed, final_psf_size);
+    psf = squeeze(sum(sum(psf, 1), 3));
     psf = psf / sum(psf(:));
 end