diff --git a/zUtil_Imaging/gtESF2LSF.m b/zUtil_Imaging/gtESF2LSF.m
index 3c5818a120471519ea4e2b5011ae8915e99ef784..0be8d4c14f468e1978b04d4ba154690d2a19a4ab 100644
--- a/zUtil_Imaging/gtESF2LSF.m
+++ b/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
diff --git a/zUtil_Imaging/gtLSFs2PSF.m b/zUtil_Imaging/gtLSFs2PSF.m
index 23d909297c257d59dd660c1403e4d640b20d591b..bd7b96ed402206d8141076ea317bae301a804240 100644
--- a/zUtil_Imaging/gtLSFs2PSF.m
+++ b/zUtil_Imaging/gtLSFs2PSF.m
@@ -1,19 +1,24 @@
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