Find stuff on an anatomical image

gofish.m

function out = gofish(mri)
% helps you find fiducial markers (realign the mri) or electrodes, assuming
% mri is a volume (3D array)
%
% use arrow keys to change slice, shift + arrow key to change +/- 10 slices
% press r to mark rpa
%       l to mark lpa
%       n to mark nasion
%       z to mark z point
%
% close figure to finish
% 
% returns struct suitable for ft_volumerealign. lpa & rpa are swapped so
% that cross(nas, zp) * rpa' > 0, allowing to resolve left & right 
% unambiguously.
%
% alternatively, use t & e to mark target and entry points for electrodes
% and a struct is returned with this information
%
% also pressing -/+ adjust upper limit on color; use shift to go faster
% look at console to see quantile used & actual clim values
%
% alternatively, press h to pop up a plot of distribution of voxel values,
% zooming and panning on this figure updates the image's contrast.
% 
% mw 13 11 2014 add pop up contrast control (type 'h')
% mw 12 11 2014 add check & resolve left-right ambiguity
% mw 02 11 2014 creation

if isstruct(mri) && isfield(mri, 'anatomy')
    anat = mri.anatomy;
else
    anat = mri;
end

f = figure(42);
f_hist = -1;
i = 1;
ijks = 'ijk';
ax = 1;
colormap(gray)
upq = 0.95;
loq = 0.05;

clim = quantile(anat(1:10:end), [0.05 upq]);
im = imagesc(squeeze(anat(i, :, :)), clim);
a = gca;

cfg = [];
cfg.method = 'fiducial';
electrodes = [];
last_electrode = '';

    % get user to click & return voxel index, taking into account which
    % axis is used for slicing
    function ijk = getijk
        figure(f)
        [ix, iy] = ginput(1);
        switch ax
            case 1, ijk = [i iy ix];
            case 2, ijk = [iy i ix];
            case 3, ijk = [iy ix i];
        end
        ijk = round(ijk);
    end

    % update slice
    function update
        figure(f)
        switch ax
            case 1, CData = anat(i, :, :); xlabel('k'), ylabel('j');
            case 2, CData = anat(:, i, :); xlabel('k'), ylabel('i');
            case 3, CData = anat(:, :, i); xlabel('j'), ylabel('i');
        end
        CData = squeeze(CData);
        XData = [1 size(CData,2)];
        YData = [1 size(CData,1)];
        set(im, 'XData', XData);
        set(im, 'YData', YData);
        xlim(XData);
        ylim(YData);
        set(im, 'CData', CData);
        title(sprintf('%s = %d', ijks(ax), i))
    end

    % reset upper quantile used for image color limit
    function req
        figure(f)
        clim = quantile(anat(1:10:end), [loq upq])
        set(a, 'CLim', clim)
    end

    function zclim(s, e)
        figure(f_hist)
        set(a, 'CLim', ylim)
        figure(f)
    end

    function key(src, ev)
        figure(f)
        % a flag used when noting landmarks
        found_point = 0;
        % shift amplifies effects
        if any(cellfun(@(c) strcmp(c, 'shift'), ev.Modifier))
            jump = 10;
        else
            jump = 1;
        end
        
        switch ev.Key
            % change slice index
            case {'uparrow', 'rightarrow'}, i = i + jump;
            case {'leftarrow', 'downarrow'}, i = i - jump;
            % change slicing dimension (i 1st, j 2nd, k 3rd)
            case {'i' 'j' 'k'}
                switch ev.Key
                    case 'i', ax = 1;
                    case 'j', ax = 2;
                    case 'k', ax = 3;
                end
                % set center slice when changing axis
                i = round(size(anat, ax) / 2);
            % close window
            case {'return', 'escape'}, close(f); return
            % change upper limit for contrast
            case {'add', 'plus'}, upq = upq * (1 + jump/100), req;
            case {'subtract', 'hyphen'}, upq = upq * (1 - jump/100), req;
            % show cumulative values allowing contrast control
            case 'h'
                if f_hist == -1
                    f_hist = figure;
                    qs = (1:1000) / 1000;
                    plot(qs, quantile(anat(1:10:end), qs));
                    set(zoom, 'ActionPostCallback', @zclim)
                    set(pan, 'ActionPostCallback', @zclim)
                    grid on
                    xlabel('quantile')
                    ylabel('voxel values')
                    title('control contrast with zoom & pan')
                    figure(f);
                else
                    figure(f_hist)
                end                
            % take note of landmarks
            case 'r', cfg.fiducial.rpa = getijk; found_point = 1;
            case 'l', cfg.fiducial.lpa = getijk; found_point = 1;
            case 'n', cfg.fiducial.nas = getijk; found_point = 1;
            case 'z', cfg.fiducial.zpoint = getijk; found_point = 1;
            % target & entry for electrodes
            case {'t' 'e'}
                ijk = getijk;
                switch ev.Key
                    case 't', lmrk = 'target'; 
                    case 'e', lmrk = 'entry'; 
                end
                prompt = ['type electrode name (' lmrk ')'];
                if ~isempty(last_electrode)
                    prompt = [prompt '(will use ''' last_electrode ''' if empty)'];
                end
                ename = input([prompt ' : '], 's');
                if isempty(ename)
                    ename = last_electrode;
                end
                last_electrode = ename;
                if strcmp(ev.Key, 't')
                    electrodes.(ename).target = ijk;
                elseif strcmp(ev.Key, 'e')
                    electrodes.(ename).entry = ijk;
                end
                found_point = 2;
            case 'shift'
                return
            otherwise
                fprintf('gofish: ignoring unknown key %s\n', ev.Key);
        end
        % update image with slice
        update;
        % print info to console if clicked for points
        if found_point == 1
            cfg.fiducial
        elseif found_point == 2
            electrodes
        end
    end

figure(f)

set(f, 'KeyPressFcn', @key);
cfg.fiducial = []
waitfor(f);

if f_hist ~= -1;
    close(f_hist)
end

if isempty(electrodes)
    if isempty(cfg.fiducial)
        return
    end
    
    if ~left_is_left(cfg.fiducial)
        warning('swapping lpa & rpa so that (nas x z+) rpa'' is positive');
        was_lpa = cfg.fiducial.lpa;
        cfg.fiducial.lpa = cfg.fiducial.rpa;
        cfg.fiducial.rpa = was_lpa;
    end
    
    out = cfg;    
else
    out = electrodes;
end

end % fidfind

function ok = left_is_left(fids)
% 
% uses fiducial points to check that left-right ambiguity is resolved 
% correctly such that (nas x z+) rpa' is positive.
%
% note: use voxel indices, not x-y-z mri coordinates!
%
% mw 12 11 2014 creation

f = fids;

% use center of pre-auricular as origin
cpa = mean([fids.rpa; fids.lpa]);
nas = fids.nas - cpa;
zp  = fids.zpoint - cpa;
rpa = fids.rpa - cpa;

%  rpa
%   |
%   |____ z+
%   \
%    \
%     nas
    
ok = cross(nas, zp) * rpa' > 0;

end % left_is_left