diff --git a/zUtil_Imaging/gtESF2LSF.m b/zUtil_Imaging/gtESF2LSF.m
new file mode 100644
index 0000000000000000000000000000000000000000..dd89ecd54b3bdbf4e9ce7ed029fe98c8a90763fa
--- /dev/null
+++ b/zUtil_Imaging/gtESF2LSF.m
@@ -0,0 +1,203 @@
+function lsf_line = gtESF2LSF(img, varargin)
+ conf = struct( ...
+ 'bb', [], 'direction', [], ...
+ 'size', 21, 'oversampling', 5, 'data_type', 'single', ...
+ 'side', 'both', 'use_astra', false, ...
+ 'keep_oversampling', true, 'verbose', true);
+ conf = parse_pv_pairs(conf, varargin);
+
+ p = gtLoadParameters();
+
+ if (ischar(img))
+ img = GtVolView.loadVolume(img);
+ elseif (numel(img) == 1)
+ 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);
+ end
+ bb(3:4) = bb(3:4) + 1 - mod(bb(3:4), 2);
+
+ conf.bb = get_roi_bb(img, 'ESF Bounding Box', bb);
+ end
+
+ if (isempty(conf.direction))
+ directions = {'left-right', 'right-left', 'top-bottom', 'bottom-top'};
+ format_dirs = {1, 2, 3, 4; directions{:}};
+ questionMsg = ['What is the rough direction of the edge (dark -> bright)?' sprintf('\n%d) %s', format_dirs{:}) '\n'];
+ dir_num = inputwdefault(questionMsg, '');
+ while (~ismember(dir_num, {'1', '2', '3', '4'}))
+ dir_num = inputwdefault(questionMsg, '');
+ end
+ conf.direction = directions{str2double(dir_num)};
+ end
+
+ if (conf.oversampling == 1)
+ edge_img = img(conf.bb(2):(conf.bb(2)+conf.bb(4)-1), conf.bb(1):(conf.bb(1)+conf.bb(3)-1));
+ else
+ shift = (1 - 1/conf.oversampling) / 2;
+ vv = linspace(conf.bb(2) - shift, conf.bb(2) + conf.bb(4) - 1 + shift, conf.bb(4) * conf.oversampling);
+ uu = linspace(conf.bb(1) - shift, conf.bb(1) + conf.bb(3) - 1 + shift, conf.bb(3) * conf.oversampling);
+ [uu, vv] = ndgrid(uu, vv);
+ edge_img = interp2(img, uu, vv, 'spline', 0);
+ end
+
+ edge_img = cast(edge_img, conf.data_type);
+
+ if (ischar(conf.direction))
+ conf.direction = find_edge_angle(edge_img, conf.direction, true);
+ end
+
+ if (conf.use_astra)
+ else
+ angle = atan2d(conf.direction(2), conf.direction(1));
+ rot_edge_img = gtRotateVolume(edge_img, -angle);
+ end
+
+ sum_esf = sum(rot_edge_img, 2)';
+ lsf = diff(sum_esf);
+ grad_pos = (1:numel(lsf)) + 0.5;
+ [max_lsf_val, max_lsf_pos] = max(lsf);
+ rel_max_lsf_pos = max_lsf_pos + 0.5;
+
+ if (conf.verbose)
+ f = figure();
+ ax = axes('parent', f);
+ plot(ax, 1:numel(sum_esf), sum_esf)
+ hold(ax, 'on')
+ plot(ax, grad_pos, lsf)
+ plot(ax, rel_max_lsf_pos, max_lsf_val, 'o')
+ hold(ax, 'off')
+ end
+
+ % Finding avg point
+ lsf_edge = (conf.size * conf.oversampling - 1) / 2;
+ inds = max_lsf_pos-5:max_lsf_pos+5;
+ center_pos = sum(lsf(inds) .* inds) / sum(lsf(inds));
+
+ lsf = interp1(lsf, linspace(center_pos-lsf_edge, center_pos+lsf_edge, conf.size * conf.oversampling), 'spline');
+ switch (conf.side)
+ case 'both'
+ case 'left'
+ lsf(lsf_edge+2:end) = flip(lsf(1:lsf_edge));
+ case 'right'
+ lsf(1:lsf_edge) = flip(lsf(lsf_edge+2:end));
+ end
+ lsf(lsf < 0) = 0;
+ if (~conf.keep_oversampling)
+ lsf = reshape(lsf, conf.oversampling, conf.size);
+ lsf = sum(lsf, 1);
+ end
+ lsf = lsf / sum(lsf);
+
+ if (conf.verbose)
+ f = figure();
+ ax = axes('parent', f);
+ plot(ax, lsf)
+ end
+
+ c = reshape(fieldnames(conf), 1, []);
+ c(2, :) = struct2cell(conf);
+ lsf_line = struct('data', lsf, c{:});
+end
+
+function bb = get_roi_bb(img, title, bb)
+ hfig = figure('name', title);
+
+ colormap(hfig, gray);
+ ax = axes('parent', hfig);
+ imagesc(img, 'parent', ax);
+ drawnow();
+ hold(ax, 'on');
+
+ plot([bb(1), bb(1)+bb(3)-1, bb(1)+bb(3)-1, bb(1), bb(1)], [bb(2), bb(2), bb(2)+bb(4)-1, bb(2)+bb(4)-1, bb(2)], 'r-');
+ bb_center = bb(1:2) + (bb(3:4) - 1) / 2;
+ plot([bb_center(1) bb_center(1)], [bb(2) bb(2)+bb(4)-1], 'c-.');
+ plot([bb(1), bb(1)+bb(3)-1], [bb_center(2) bb_center(2)], 'c-.');
+
+ % interactive check
+ questionMsg = 'Are you satisfied with this bounding box for the region where to compute the ESF sampling? [y/n]';
+ while (strcmpi(inputwdefault(questionMsg, 'y'), 'n'))
+ disp('Please select (zoom onto) the region where to compute the ESF and disable zoom afterwards');
+
+ figure(hfig), clf;
+ ax = axes('parent', hfig);
+ imagesc(img, 'parent', ax);
+ h = zoom(hfig);
+ set(h, 'Enable', 'on');
+ waitfor(h, 'Enable', 'off')
+ hold(ax, 'on');
+ disp('-> Now click top-left and bottom-right corners for final selection');
+ bb = ginput(2);
+ bb = round([bb(1, 1), bb(1, 2), bb(2, 1)-bb(1, 1)+1, bb(2, 2)-bb(1, 2)+1]);
+
+ % Make sure the bbox width and height are odd numbers (needed for correlation)
+ bb(3:4) = bb(3:4) + 1 - mod(bb(3:4), 2);
+
+ % Plot the new bounding box
+ plot([bb(1), bb(1)+bb(3)-1, bb(1)+bb(3)-1, bb(1), bb(1)], [bb(2), bb(2), bb(2)+bb(4)-1, bb(2)+bb(4)-1, bb(2)], 'r-');
+ bb_center = bb(1:2) + (bb(3:4) - 1) / 2;
+ plot([bb_center(1) bb_center(1)], [bb(2) bb(2)+bb(4)-1], 'c-.');
+ plot([bb(1), bb(1)+bb(3)-1], [bb_center(2) bb_center(2)], 'c-.');
+ end
+
+ close(hfig);
+
+ fprintf('Chosen BB: [%s]\n', sprintf(' %d', bb))
+end
+
+function [direction, angle] = find_edge_angle(edge_img, rough_dir, do_plot)
+ % We rotate the roi, to facilitate the computation
+ switch (rough_dir)
+ case 'left-right'
+ base_k_90_rot = 0;
+ case 'right-left'
+ base_k_90_rot = 2;
+ case 'top-bottom'
+ base_k_90_rot = 3;
+ case 'bottom-top'
+ base_k_90_rot = 1;
+ otherwise
+ error('gtESF2LSF:wrong_argument', ...
+ 'Direction "%s" unknown', rough_dir);
+ end
+ edge_img = rot90(edge_img, base_k_90_rot);
+
+ % We will now use the derivative of all the lines, to determine the
+ % angle
+ grads = diff(edge_img, 1, 1);
+ [~, max_grads_pos] = max(grads, [], 1);
+
+ % least square fitting of the tilt angle
+ X = [ones(numel(max_grads_pos), 1), (1:numel(max_grads_pos))'];
+ vals = X \ reshape(max_grads_pos, [], 1);
+
+ if (exist('do_plot', 'var') && do_plot)
+ f = figure();
+ ax = axes('parent', f);
+ plot(ax, max_grads_pos)
+ hold(ax, 'on')
+ plot(ax, (1:numel(max_grads_pos))', X*vals);
+ hold(ax, 'off')
+ end
+
+ angle = base_k_90_rot * 90 + atand(vals(2));
+
+ direction = [cosd(angle), sind(angle)];
+end
diff --git a/zUtil_Imaging/gtLSFs2PSF.m b/zUtil_Imaging/gtLSFs2PSF.m
new file mode 100644
index 0000000000000000000000000000000000000000..23d909297c257d59dd660c1403e4d640b20d591b
--- /dev/null
+++ b/zUtil_Imaging/gtLSFs2PSF.m
@@ -0,0 +1,95 @@
+function psf = gtLSFs2PSF(lsfs, varargin)
+ conf = struct('method', 'fourier');
+ conf = parse_pv_pairs(conf, varargin);
+
+ num_lsfs = numel(lsfs);
+
+ 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);
+ 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 = 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;
+ case 'fourier'
+ data_f = ifftshift(data);
+ data_f = fft(data_f);
+ data_f = fftshift(data_f);
+
+ psf_edge = (numel(lsfs.data) - 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);
+
+ psf_f = ifftshift(psf_f);
+ psf = ifft2(psf_f);
+ psf = real(fftshift(psf));
+ case 'hankel'
+ end
+ else
+ % Checking compatibility
+ if (any(lsfs(1).size ~= [lsfs(2:end).size]))
+ error('gtPSFEstimateFromRef:wrong_argument', ...
+ 'All LSFs should have the same size');
+ end
+ if (any([lsfs.keep_oversampling]))
+ if (any(lsfs(1).keep_oversampling ~= [lsfs(2:end).keep_oversampling]))
+ error('gtPSFEstimateFromRef:wrong_argument', ...
+ 'All LSFs should behave the same regarding oversampling');
+ end
+ if (any(lsfs(1).oversampling ~= [lsfs(2:end).oversampling]))
+ error('gtPSFEstimateFromRef:wrong_argument', ...
+ 'All LSFs should have the same oversampling');
+ end
+ 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;
+
+ dirs = cat(1, lsfs.direction);
+ angles = mod(atan2d(dirs(:, 2), dirs(:, 1)), 360);
+
+ % Let's bring everything in a polar re-arrangement
+ data_f = [data_f(:, 1+psf_edge:end); flip(data_f(:, 1:1+psf_edge), 2)];
+ angles = [angles; mod(angles - 180, 360)];
+
+ % let's re-order them
+ [angles, inds] = sort(angles);
+ data_f = data_f(inds, :);
+
+ % Let's wrap around: we need one going to more than 180, and one
+ % going to less than -180, and left and right sides have to be
+ % treated as radial
+ angles = [angles(end) - 360; angles; angles(1) + 360];
+ data_f = [data_f(end, :); data_f; data_f(1, :)];
+
+ [xx, yy] = ndgrid(-psf_edge:psf_edge, -psf_edge:psf_edge);
+ rr = sqrt(xx .^ 2 + yy .^ 2);
+ aa = mod(atan2d(yy, xx), 360);
+
+ [aa_base, rr_base] = ndgrid(angles, 0:psf_edge);
+ psf_f = interp2(rr_base, aa_base, data_f, rr, aa, 'spline', 0);
+
+ psf_f = ifftshift(psf_f);
+ psf = ifft2(psf_f);
+ 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));
+ end
+ psf = psf / sum(psf(:));
+end