function  [layout,locate] = biac_montage(Y,plane,FS,CloseUp)
% FORMAT  [ layout, locate ] = biac_montage(Y,plane,gap,FS,CloseUp)
%   Called by function biac_movie.
%
% FUNCTION:
%   Convert the 3D scan image Y into a 2D montage image of slices.
%   Also returns vector "locate" to invert the process: from pixel -> voxel.
%
% INPUTS:  
%   Y is 3D matrix for a single image, comes in as type double.
%   plane is the plane to display
%   FS is the first slice to display in the CloseUp=1 mode
%   CloseUp = 0 gives all slices, with one pixel/voxel
%          = 1 gives 20-25 consecutive slices, with four pixels/voxel.
% OUTPUTS: 
%  layout is a 2D image of all slices in designated plane, as type double.
%  locate is a length 6 vector used to locate the voxel from a pixel in layout.
%
% NOTE:  Large image arrays will be clipped to size (80,100,80).
%
%See also:   display_slices, show_image, slice_overlay,  montage (in Matlab toolbox)

%Note: There may be a faster way to do this using permute and reindexing.

% Maximum sizes allowed. Function will clip number and size of slices to fit.
[ A ] = size(Y);
NX = A(1);
NY = A(2);
NZ = A(3);
NX = min(NX,80);   % Slices are clipped for images larger than 80x100x80.
NY = min(NY,100);  %   Bug: Would be nicer to clip in the middle.
NZ = min(NZ,80);
gap = 0;           %  No gaps allowed between slices.
%  Y(5,10,15) = 2000;   % Test Trace point for unmontage function.

% Make an image array
if CloseUp == 1  % For the 20-25 slice montage case
    % Three cases for slice sizes of the large image (79,95,69)
    %   Coronal   25 slices of size 79x70.  Montage is  5x5, size 790x700. 
    %   Sagittal  24 slices of size 95x70.  Montage is  4x6, size 760x840.
    %   Axial     20 slices of size 95x79.  Montage is  4x5, size 790x760
    if plane == 1  % axial/transverse    % was 4,5
        jlim = NX;  klim = NY;  NROWS = 5;  NCOLS = 4;  LastSlice = min(NZ,FS+19);
    elseif  plane == 2  %  sagittal
        jlim = NY;  klim = NZ;  NROWS = 4;  NCOLS = 6;  LastSlice = min(NX,FS+23);
    else   %  plane == 3  % coronal
        jlim = NX;  klim = NZ;  NROWS = 5;  NCOLS = 5;  LastSlice = min(NY,FS+24);
    end 
    layout = zeros(2*NROWS*jlim+6*gap,2*NCOLS*klim+6*gap);  % for the montage
    temp = zeros(jlim,klim);  % for a single slice
else %  CloseUp = 0   %  For viewing all the slices
    % Three cases for large image (79,95,69)
    %   Coronal   95(NY) slices of size 79x70.  Montage is 10x10, size 790x700. 
    %   Sagittal  79(NX) slices of size 95x70.  Montage is  8x10, size 760x700.
    %   Axial     69(NZ) slices of size 95x79.  Montage is  7x10, size 665x790
    if plane == 1  % axial/transverse
        jlim = NX;  klim = NY;  NROWS = 10;  NCOLS = 7;  LastSlice = NZ;  % was 7,10
    elseif  plane == 2  %  sagittal
        jlim = NY;  klim = NZ;  NROWS = 8;  NCOLS = 10;  LastSlice = NX;
    else   %  plane == 3  % coronal
        jlim = NX;  klim = NZ;  NROWS = 10; NCOLS = 10;  LastSlice = NY;
    end
    layout = zeros(NROWS*jlim,NCOLS*klim);  
    temp = zeros(jlim,klim);
    FS = 1;
end

if CloseUp == 0    % All slices
	for icolumn = 1:NCOLS
        fsn = FS + NROWS*(icolumn-1);  % first slice to display
		for i = fsn:fsn + NROWS - 1
            if i <= LastSlice   % Montage is larger than number of slices
                if (plane == 2) temp(:,:) = Y(NX-i+1,1:jlim,1:klim); end
                if (plane == 3) temp(:,:) = Y(1:jlim,NY-i+1,1:klim); end
                if (plane == 1) temp(:,:) = Y(1:jlim,1:klim,NZ-i+1); end
                ix = gap + (i - fsn )*(jlim+gap);   % depend on (i-fsn)
                iy = gap + (icolumn -1)*(klim+gap);    % depend on icolumn
                for j = 1:jlim
                    for k = 1:klim   % for all slices,  use one pixel per voxel
                        layout(ix+j  ,iy+k)   = temp(j,k);
                    end
                end
            end
        end
    end
end

if CloseUp == 1  %  CloseUp version
	for icolumn = 1:NCOLS
        fsn = FS + NROWS*(icolumn-1);  % first slice to display
		for i = fsn:fsn + NROWS - 1
            if i <= LastSlice   % Requested slices must be in the data set.
                if (plane == 2) temp(:,:) = Y(NX-i+1,1:jlim,1:klim); end
                if (plane == 3) temp(:,:) = Y(1:jlim,NY-i+1,1:klim); end
                if (plane == 1) temp(:,:) = Y(1:jlim,1:klim,NZ-i+1); end
                ix = gap + (i - fsn )*(2*jlim+gap);   % depend on (i-fsn)
                iy = gap + (icolumn -1)*(2*klim+gap);    % depend on icolumn
                for j = 1:jlim
                    for k = 1:klim   % for CloseUp, use four pixels per voxel
                        layout(ix+2*j  ,iy+2*k)   = temp(j,k);
                        layout(ix+2*j+1,iy+2*k)   = temp(j,k);
                        layout(ix+2*j  ,iy+2*k+1) = temp(j,k);
                        layout(ix+2*j+1,iy+2*k+1) = temp(j,k);
                    end
                end
            end
		end
    end
end

% Arrange the image to get these orientations:
%    Axial - eyes up in a slice, slices with eyes in upper images
%    Sagittal - brain stem on right side
%    Coronal - slices with eyes in lower images.
layout = layout';  
layout = flipud(layout);
layout = fliplr(layout);

%  Return values used to locate a voxel from a pixel in the function biac_unmontage.
locate(1) = plane;
locate(2) = NCOLS; %  nxrect, number slices in i-direction
locate(3) = NROWS; %  nyrect, number slices in j-direction
locate(4) = klim;  %  xrect, slice size in i-direction
if (CloseUp == 1) locate(4) = 2*locate(4); end
locate(5) = jlim;  %  yrect, slice size in j-direction
if (CloseUp == 1 ) locate(5) = 2*locate(5); end
locate(6) = LastSlice;   %

% Return to main function with layout array and locate parameters.

% NOTES: Procedure to return an (X,Y,Z) voxel from a pixel location.
%   Implemented in function "unmontage".
% For the calculation in this function:
%   plane = 3 coronal, slice size is NZxNX, arrayrect(10,10)
%   plane = 2 sag., slice size is NZxNY, arrayrect(10,8)
%   plane = 1 axial, slice size is NYxNX, arrayrect(7,10)
%   If CloseUp = 1, need to double jlim and klim, and use LastSlice.
% To invert the process, need the following steps:
%   For the plane orientation, need to determine slice #, x and y values.
%   xrect, yrect is slice size in (i,j) direction.
%   nxrect, nyrect is number of slices in (i,j) direction
% For a pixel (I,J) is layout:
%   xc = floor(I/locate(4)) + 1;
%   yc = floor(j/locate(5)) + 1;
%   snum = xc*locate(3) - (locate(2) - yc)
%   if ( snum < 1 ) snum = 1; end
%   jc = J mod(locate(5));   %yrect
%   ic = I mod(locate(4)); % xrect
%   jca = yrect - jc;
%   ica = xrect - ic;
%   if plane = 3,  X = jca,  Z = ica, Y = slice
%   if plane = 2,  Y = jca,  Z = ica, X = slice
%   if plane = 1,  X = jca,  Y = ica, Z = slice

