% SHOW_STACK   Show maximum intensity projection of an image stack.
% (trees package)
%
% HP = show_stack (stack, options)
% --------------------------------
%
% Show the maximum intensity projections of a stack on a 3D patch (see
% "load_stack" for a description of a stack structure and see "cgui_tree"
% stk_ panel).
%
% Inputs
% ------
% - stack::stack structure: (see load_stack)
% - options::string: {DEFAULT: '-a'}
%     '-a' : axis equal
%
% Output
% ------
% - HP::handle: graphical handle to the graphics objects
%
% Example
% -------
% HP = show_stack (stack);
%
% See also load_stack
% Uses
%
% the TREES toolbox: edit, visualize and analyze neuronal trees
% Copyright (C) 2009  Hermann Cuntz

function HP = show_stack (stack, options)

if (nargin<2)||isempty(options),
    options = '-a'; % {DEFAULT: axis equal}
end

HP = zeros (1, length (stack.M) * 3);
for ward = 1 : length (stack.M),
    HP ((ward - 1) * 3 + 1) = surface ( ...
        -.5 * stack.voxel (1) + (stack.coord (ward, 1) + ...
        stack.voxel (1) * ...
        [0 size(stack.M {ward}, 2)-1; 0 size(stack.M {ward}, 2)-1]), ...
        -.5 * stack.voxel (2) + (stack.coord (ward, 2) + ...
        stack.voxel (2) * ...
        [size(stack.M {ward}, 1)-1 size(stack.M {ward}, 1)-1; 0 0]), ...
        (stack.coord (ward, 3) + repmat (0, 2, 2)));
    set (HP ((ward - 1) * 3 + 1), ...
        'CData', flipud (double (max (stack.M {ward}, [], 3))), ...
        'FaceColor', 'texturemap', 'Edgecolor', 'none', 'facealpha', 0.5);
    HP ((ward - 1) * 3 + 2) = surface ( ...
        -.5 * stack.voxel (1) + (stack.coord (ward, 1) + ...
        stack.voxel (1) * ...
        [0 size(stack.M {ward}, 2)-1; 0  size(stack.M {ward}, 2)-1]), ...
        (stack.coord (ward, 2) + repmat (0, 2, 2)), ...
        -.5 * stack.voxel (3) + (stack.coord (ward, 3) + ...
        stack.voxel (3) * ...
        [size(stack.M {ward}, 3)-1 size(stack.M {ward}, 3)-1; 0 0]));
    set (HP ((ward - 1) * 3 + 2), ...
        'CData', flipud (double (squeeze (max (stack.M {ward}, [], 1)))'), ...
        'FaceColor', 'texturemap', 'Edgecolor', 'none','facealpha',0.5);
    HP ((ward - 1) * 3 + 3) = surface ( ...
        (stack.coord (ward, 1) + repmat (0, 2, 2)), ...
        -.5 * stack.voxel (2) + (stack.coord (ward, 2) + ...
        stack.voxel (2) * ...
        [0 size(stack.M {ward}, 1)-1; 0  size(stack.M {ward}, 1)-1]), ...
        -.5 * stack.voxel (3) + (stack.coord (ward, 3) + ...
        stack.voxel (3) * ...
        [size(stack.M {ward}, 3)-1 size(stack.M {ward}, 3)-1; 0 0]));
    set (HP ((ward - 1) * 3 + 3), ...
        'CData', flipud (double (squeeze (max (stack.M {ward}, [], 2)))'), ...
        'FaceColor', 'texturemap', 'Edgecolor', 'none', 'facealpha', 0.5);
end

if strfind (options, '-a'),
    axis equal;
end