PSF: Added initial functions for dealing with ESF, LSF and producing a PSF.

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

zUtil_Imaging/gtESF2LSF.m

+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

zUtil_Imaging/gtLSFs2PSF.m

+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