function h = slice3i(vol, I2X, slicedim, sliceidx, handle, keepxyz)
% Display a slice from a volume in 3-D
% h = slice3(vol, I2X, slicedim, sliceidx, handle)
%
% Vol is either a scalar or RGB volume, e.g. N x M x K or N x M x K x 3.
% I2X is a transformation matrix from volume coordinates to xyz 3-D world
% coordinates, similar to the transform used by image3i.
%
%     [x y z 1]' = I2X * [i j k 1]'
%
% slicedim  is 1, 2 or 3 and corresponds to i, j and k.
% sliceidx  the index of the slice along slicedim
% handle    an optional handle to a previous slice to reuse
% h         the handle to the created slice.
%
% Example:
%
% load mri;
% T = [1 0 0 0;0 1 0 0;0 0 2.5 0];
% h1 = slice3i(squeeze(D),T,1,64);
% h2 = slice3i(squeeze(D),T,2,64);
% h3 = slice3i(squeeze(D),T,3,14);
% colormap gray(88);
% view(30,45); axis equal;
%
% SEE ALSO: image3i, imagesc3
%
% Author: Anders Brun, anders@cb.uu.se (2009)
%
% Modified by Qianqian Fang, q.fang at neu.edu
%
% Features added:
%    1) compatible with octave 4.0 or newer
%    2) start 3D rotation by clicking/holding the middle mouse key
%    3) adjust colormap levels by dragging with the right mouse key
%    4) support 4-D data (up/down key to change frames)
%

if (nargin < 6)
    keepxyz = 0;
end

try
    h = update_slice(vol, I2X, slicedim, sliceidx, handle, keepxyz);
catch
    h = update_slice(vol, I2X, slicedim, sliceidx);
end

% set up gui
gui.handle = h;
gui.vol = vol;
gui.I2X = I2X;
gui.slicedim = slicedim;
gui.sliceidx = sliceidx;
set(gui.handle, 'ButtonDownFcn', @startmovit);

set(h, 'UserData', gui);

function stoprotate(src, hrot)
if strcmp(get(gcf, 'SelectionType'), 'extend')
    rotate3d off;
end

function startmovit(src, evnt)

if strcmp(get(gcf, 'SelectionType'), 'extend')
    rotate3d on;
    return
end

% Unpack gui object
gui = get(src, 'UserData');
% turn off mouse pointer
% set(gcf,'PointerShapeCData',nan(16,16));
% set(gcf,'Pointer','custom');

thisfig = gcbf();
gui.startray = get(gca, 'CurrentPoint');
gui.startidx = gui.sliceidx;
gui.p0 = get(gcf, 'currentpoint');
gui.level0 = size(colormap, 1);

% Store gui object
set(src, 'UserData', gui);
set(thisfig, 'WindowButtonMotionFcn', @movit);
set(thisfig, 'WindowButtonUpFcn', @stopmovit);
set(thisfig, 'UserData', src);

function movit(src, evnt)
% Unpack gui object
gui = get(get(gcf, 'UserData'), 'UserData');
% Some safetymeasures
try
    if isequal(gui.startray, [])
        return
    end
catch
end

if strcmp(get(gcf, 'SelectionType'), 'alt')
    if (isempty(gui) || ~isfield(gui, 'p0'))
        return
    end
    delta = get(gcf, 'currentpoint') - gui.p0;
    windim = get(gcf, 'position');
    delta = delta ./ windim(3:4);
    currentlevel = round(gui.level0 + delta(2) / 4 * 256);
    if (currentlevel > 256)
        currentlevel = 256;
    elseif (currentlevel < 8)
        currentlevel = 8;
    end
    colormap(jet(currentlevel));
    return
end

% Do "smart" positioning of the markers...
nowray = get(gca, 'CurrentPoint');

% Project rays on slice-axis
s = gui.I2X(1:3, gui.slicedim);
a = gui.startray(1, :)';
b = gui.startray(2, :)';

alphabeta = pinv([s' * s, -s' * (b - a); (b - a)' * s, -(b - a)' * (b - a)]) * [s' * a, (b - a)' * a]';
pstart = alphabeta(1) * s;
alphastart = alphabeta(1);
a = nowray(1, :)';
b = nowray(2, :)';
alphabeta = pinv([s' * s, -s' * (b - a); (b - a)' * s, -(b - a)' * (b - a)]) * [s' * a, (b - a)' * a]';
pnow = alphabeta(1) * s;
alphanow = alphabeta(1);
slicediff = alphanow - alphastart;

gui.sliceidx = gui.startidx + slicediff;
gui.sliceidx = min(max(1, gui.sliceidx), size(gui.vol, gui.slicedim));
update_slice(gui.vol, gui.I2X, gui.slicedim, gui.sliceidx, gui.handle);
drawnow;

% Store gui object
set(get(gcf, 'UserData'), 'UserData', gui);

function stopmovit(src, evnt)

thisfig = gcbf();
% set(gcf,'Pointer','arrow');
set(thisfig, 'WindowButtonUpFcn', '');
set(thisfig, 'WindowButtonMotionFcn', '');
gui = get(get(gcf, 'UserData'), 'UserData');
gui.startray = [];
% gui.sliceidx = gui.nextsliceidx
set(get(gcf, 'UserData'), 'UserData', gui);
set(gcf, 'UserData', []);
drawnow;

function h = update_slice(vol, I2X, slicedim, sliceidx, handle, keepxyz)

if (nargin < 6)
    keepxyz = 0;
end

if ndims(vol) == 3         % Scalar mode
elseif ndims(vol) == 4     % RGB mode
    if (size(vol, 4) > 3)
        gui = get(get(gcf, 'UserData'), 'UserData');
        if (~isempty(gui) && isfield(gui, 'frame'))
            vol = gui.vol(:, :, :, gui.frame);
        else
            vol = gui.vol(:, :, :, 1);
        end
    end
else
    error('Only scalar and RGB images supported');
end

% Create the slice
if slicedim == 3 % k
    ij2xyz = I2X(:, [1 2]);
    ij2xyz(:, 3) = I2X * [0 0 sliceidx 1]';
    sliceim = squeeze(vol(:, :, round(sliceidx), :));

elseif slicedim == 2 % j
    ij2xyz = I2X(:, [1 3]);
    ij2xyz(:, 3) = I2X * [0 sliceidx 0 1]';
    sliceim = squeeze(vol(:, round(sliceidx), :, :));

elseif slicedim == 1 % i
    ij2xyz = I2X(:, [2 3]);
    ij2xyz(:, 3) = I2X * [sliceidx 0 0 1]';
    sliceim = squeeze(vol(round(sliceidx), :, :, :));
else
    error('Slicedim should be 1, 2 or 3');
end

if nargin < 5 || handle == 0
    h = image3i(sliceim, ij2xyz);
else
    h = image3i(sliceim, ij2xyz, handle, keepxyz);
end