function biac_movie(Action)
% FORMAT biac_movie
%      A command line user function. All inputs are asked by GUI.
%
% FUNCTIONS:  Displays bulk data for quick visual review by a user. Typical
%   time to process and display all the data in 100 scans of size (79,95,69)
%   is a few minutes. Every voxel of scan data can be seen.
%       Each scan is made into a montage of slices with orientation
%   chosen by input. The montage may be all slices or about 20 consecutive
%   slices. The slice displays may either be image data or contrast data.
%   Contrast data shows the difference between each image and a
%   reference image, amplified so that small data variations are more visible.
%       In movie mode, every scan montage is shown, followed by a time history
%   of mean intensity and position of each scan. In slider mode, every scan
%   montage is available by selecting the slider position. The display is
%   a default Matlab window where the zoom button can be used
%   to examine individual pixels. 
%    
% INPUTS
%   Choose Orientation:  Transverse, Sagittal, or Coronal.
%   Choose a Set of Images (in .img format) to review.
%   Choose a Range of images to process. Default limit is 100 total scans.
%   Select All slices, or Some slices.
%       If some slices, select the first slice to display. About 20-25
%       consecutive slices are shown, at twice the image size per slice.
%   Select Raw or Contrast display. 
%       Raw shows the image data in the set of scans. 
%       Contrast shows the difference of the image from a reference scan. The
%          reference is chosen 1/3 of the way through the range of scans.
%          The display shows yellow as higher, and blue as lower from the 
%          reference. Values are amplified to be 5x more sensitive for viewing.
%          Contrast mode helps to see small errors in processed images.
%   Select Movie & Time History,  or interactive Slider.
%       If a movie, select a movie frame rate. Three loops are played.
%          1 or 2 fps is slow enough to visually spot artifacts.
%          At the end of the movie, a time history of the average intensity 
%          and position of each scan is displayed.
%       The slider shows individual frames and allows zoom in and out.
%           
% OUTPUTS
%   The reference image, automatically chosen 1/3 of the way through the selected range.
%   A movie of a montage of slices, for all the scans of data, 
%          followed by a time history plot, OR
%   An interactive slider display with all frames of the montage movie.
%
% BUGS
%    Left-Right is not well-defined on the output images.
%    Out of Memory error message may occur for large numbers of scans.
%        Reduce the number of scans, or close other running programs.
%    The reference scan will show as solid black in contrast viewing.
%
% See also:  spm_movie,  AnalyzeMovie
%
% Paul Mazaika  April 2004.  Derived as a faster version of spm_movie.

% SIZE LIMIT - This parameter can be adjusted for a particular computer.
%   One volume Y of size (79,95,68) is 4 MB when read as type double.
%   100 volumes of Y is 50 MB when stored at type uint8.
%   An experiment session may contain 300-600 volumes.
MAXSCAN = 100;

if (nargin==0)  % For calling in line with arguments
    Orient= spm_input('Select slice orientation',...
		1,'m','transverse|sagittal|coronal');
	F      = spm_get(Inf,'.img','select images');
    scansel = spm_input('Select range of up to 100 scans e.g. 11:110',1);  %array of integers
    %CloseUp = spm_input('Select All Slices(0) or 25 slices(1)',1);
    CloseUp = spm_input('All Slices or 25 larger slices',1, ' All | 25 larger', [ 0 1 ] );
     if CloseUp == 1
         slicesel = spm_input('Select first slice of 25 slices, e.g. 11',1);   % array of integers
     else
         slicesel = 1;    % Start at first slice for All slice case
     end
    typesel = spm_input('Select Raw or Contrast data',1,'b',' Raw | Contrast', [ 1 2 ] );
    smode=spm_input('Select viewing mode',1,'b',' Movie | Slider', [ 0 1 ] );
    %smode=spm_input('View movie(m) or Use Slider(s) ? (m/s) ',1,'s','');
        if smode == 1   % slider
           nloops=1;
           mode=1;
        else   %  smode == 0, movie
           frate = spm_input('Frames/second ? (e.g. 2)',1);
           nloops = 3;   %  spm_input('# loops',2);
           mode=0;
        end
%  Preliminary sizing
     f = deblank(F(1,:));
     V = spm_vol(f);
     nFsize = size(F,1);
     Action = 'Load'; nnnn = 6;
end
switch lower(Action), case('load')
    
%  ---------- Load the Data -----------------------------------------
% Last scan is minimum of user selection, size limited, or end of group.
% Note Scan 1 is first scan in group; its name is usually not image_001. !!!
firstscan = scansel(1);   
lastscanU = firstscan + length(scansel) - 1;  % User choice for last scan
lastscanM = firstscan + MAXSCAN - 1;  % Size limit for memory and program speed
lastscanF = nFsize;       % Size limit by end of available files.
lastscan = min( [ lastscanU lastscanM lastscanF ] );
Slimit = lastscan - firstscan + 1;   % Number of scans that will be processed.

% Choose preview image about 1/3 of the way in, far from initialization.
nbase = round(firstscan + (lastscan-firstscan)/3);
f = deblank(F(nbase,:));
V = spm_vol(f);
gap = 0;   %  No space allowed between image slices in the montage.
% Read and draw the preview image.
Y = spm_read_vols(V);   % Y is type double, size 4MB for (79,95,68); 
% Y is a 4D representation of images, 4th dimension is the scan number.
[layout,nil] = biac_montage(Y,Orient,slicesel,CloseUp);    % Currently double
figure(4);   % Standard Matlab window has zoom enlargement, etc.
imagesc(layout,[ 50 2000 ]); colormap(gray)   % puts one image up to start.
drawnow;
disp('Reference image has been read and displayed.')

% Set up the arrays. Size of the image was determined in function biac_montage.
[ xs, ys ] = size(layout);
    % In movie mode, data is stored in movie frames, not in the imageblock array.
        if (mode == 0)  Slimit = 1; end
    % Dimension imageblock on the fly, else "zeros" command makes type double.
    % imageblock = uint8(zeros(xs,ys,Slimit));  % Size is 50MB for 100 scans of uint8 data.
    temp = uint16(zeros(xs,ys));
    temp1 = zeros(xs,ys);
    layout8 = uint8(zeros(xs,ys)); 
    history = zeros(lastscan,4);  %  Track some global properties 
    
% SET THE COLORMAP AND CONTRAST 
% For a raw image, Range 0-2040 is mapped to range of 64.
% For contrast image, amplitude scale will be 5x larger. -160,+160 -> 0,64.   
cmap = colormap(gray);   % Size cmap is 64, so must map imageblock values to 0-63.
if ( typesel == 2 )   % Get the contrast map.
    cmap = biac_blue2yellowmap(1);   % cmap = colormap(cool) is also pretty good.
end

% Load the data, flatten into a montage, and convert into uint8.
% Reads a scan at a time, and sets up a uint8 image, or a movie frame, for it.
disp('Loading data - this may take a few minutes for a lot of scans')
for iscan = firstscan:lastscan
    nscan = iscan - firstscan + 1;
    f =deblank(F(iscan,:));
    V = spm_vol(f);
    Y = spm_read_vols(V);     
    history(iscan,1:4) = biac_GlobalProp(Y);  %  Collect time histories
    [temp1,nil] = biac_montage(Y,Orient,slicesel,CloseUp);  % temp is double 
    if mode == 1   %  Slider mode-  load imageblock
        if typesel == 1   % For raw data, we'll divide by 64 for the movie, anyway.
            temp = uint16(temp1);   % Prepares non-integer data for bitshift operation.
            imageblock(:,:,nscan) = uint8( bitshift(uint8(bitshift(temp,-3)),-2)) ;   % Range 0-2040
        else  % typesel == 2     % For contrast data, compare to preview scan.
            temp = uint16(max(1.6*(temp1-layout)+254,0));  % clips the bottom
            temp = min(temp,510);  % clips the top
            imageblock(:,:,nscan) = uint8( bitshift(temp,-3));   % Range -250 to +250 -> 64.
        end
    else   %   mode == 0,   Movie mode- load the movie image frames
        if typesel == 1   % For raw data, we'll divide by 64 for the movie, anyway.
            temp = uint16(temp1);   % Prepares non-integer data for bitshift operation.
            imageblock(:,:,1) = uint8( bitshift(uint8(bitshift(temp,-3)),-2)) ;   % Range 0-2040
        else  % typesel == 2     % For contrast data, compare to reference scan.
            temp = uint16(max(1.6*(temp1-layout)+254,0));  % clips the bottom
            temp = min(temp,510);  % clips the top
            imageblock(:,:,1) = uint8( bitshift(temp,-3));   % Range -250 to +250 -> 64.
        end
        layout8 = min(63,imageblock);  %  Guarantees it fits into movie uint8 format.
        M(nscan) = im2frame(layout8,cmap);
    end
end

% If there is a lot of storage, it might be faster to read the entire multi-scan image structure
%   at once, and then create the images. But storage req'ts will be doubled.


% Display movie or slider as required
%----------------------------------------------------------
ha =figure('Position',[100 100 830 830]);  % Create the movie or slider figure.
figure(ha);
colormap(cmap);
if mode==0   %  Movie Mode and Time History
   h=findobj('Tag','ToDelete');
   if ~isempty(h)
      delete(h);
   end;
   movie(ha,M,nloops,frate,[15 15 0 0])   % Default is frate = 4 frames/sec.
   %  Show the time history plots at the end of the movie.
   subplot(4,1,1);
   plot(history(firstscan:lastscan,1));
   title('Scan History');
   ylabel('Intensity');
   subplot(4,1,2);
   plot(history(firstscan:lastscan,2)); 
   limlow = mean(history(firstscan:lastscan,2)) - 0.3;  % 30% voxel from the mean
   ylim([ limlow  limlow+0.6 ]);
   ylabel('x-movement (voxels)');
   subplot(4,1,3);
   plot(history(firstscan:lastscan,3));
   limlow = mean(history(firstscan:lastscan,3)) - 0.3;
   ylim([ limlow  limlow+0.6 ]);
   ylabel('y-movement (voxels)');
   subplot(4,1,4);
   plot(history(firstscan:lastscan,4));
   limlow = mean(history(firstscan:lastscan,4)) - 0.3;
   ylim([ limlow  limlow+0.6 ]);
   ylabel('z-movement (voxels)');
   xlabel(['Relative Scan number. Relative Scan 1 is Scan ' num2str(firstscan-1)]);
elseif mode==1   %  Slider Mode
   if Slimit == 1   %  Slimit is number of scans processed.
      min_step=0.5;
      max_step=1;
      min_slide=0;
      max_slide=0.5;
   elseif Slimit == 2
      min_step=0.5;
      max_step=1;
      min_slide=0;
      max_slide = Slimit - 1;
   elseif Slimit < 10
      min_step=1/(Slimit-1);
      max_step=1;
      min_slide=0;
      max_slide = Slimit - 1;
   else
      min_step=1/(Slimit-1);
      max_step=10/(Slimit-1);
      min_slide = 0;
      max_slide = Slimit - 1;
   end
   h=findobj('Tag','ToDelete');
   if ~isempty(h)
      delete(h);
   end;
   figwin = ha;
   m_struct=struct('movie',imageblock,'filename',F,'first',firstscan,...
       'typesel',typesel,'orient',Orient);
   s=uicontrol('Style','slider','Parent',figwin,...      
		    'Position',[245 80 338 30],...
		    'Min',min_slide,'Max',max_slide,...
		    'SliderStep',[min_step max_step],...
                    'Callback','biac_movie(''Scroll'');',...
		    'Userdata',m_struct,'String','Frame number');
   if typesel == 2   % Draw the colormap for the contrast case.
       figbar = axes('position',[ 0.93 0.3 0.03 0.4]);
       darray = [63:-1:0];
       image(darray');
       set(figbar,'Ytick',[1 17 33 48 64],'YTickLabel',[+160 +80 0 -80 -160],'XTickLabel',[]);
   end
   set(0,'CurrentFigure',figwin);
   set(0,'ShowHiddenHandles','on')%%%%%%%%%%%%I added this
   fig = axes('position',[0.07 0.15 0.8 0.8],'Parent',figwin);
   I = imageblock(:,:,1);
   image(I,'Parent',fig,'Tag','ToDelete');
   axis image;
   axis off;
   frame=sprintf('%s',deblank(F(firstscan,:)));  % Gets the path and file name
   t=title(frame,'FontSize',14,'Interpreter','none','Tag','ToDelete');
end
case('scroll')
%==========================================================
global PRINTSTR
if isempty(PRINTSTR)
   PRINTSTR='print -dpsc2 -append spm.ps';
end;
g=get(gcbo,'Value');
m_struct=get(gcbo,'Userdata');
imageblock=m_struct.movie;
Or=m_struct.orient;
firstscan=m_struct.first;
typesel = m_struct.typesel;
F=m_struct.filename;
[ htemp, currfig ] = gcbo;
figwin= figure(currfig);  % spm_figure('FindWin','Graphics');
if typesel == 2   % Draw the colorbar for the contrast case.
    figbar = axes('position',[ 0.93 0.3 0.03 0.4]);
    darray = [63:-1:0];
    image(darray');
    set(figbar,'Ytick',[1 17 33 48 64],'YTickLabel',[+160 +80 0 -80 -160],'XTickLabel',[]);
end
fig = axes('position',[0.07 0.15 0.8 0.8],'Parent',figwin);
% if isempty(figwin)
%    figwin=spm_figure('Create','Graphics','Graphics');
% end
h=findobj('Tag','ToDelete');
%%H  = gca;
%j=get(H,'children');
if ~isempty(h)
    delete(h);
    %delete(j(1,1))
end;

figure(currfig);
nscan = size(imageblock,3);
gscan = 1 + floor(g+0.5);
if (gscan > nscan) gscan = nscan; end
I=imageblock(:,:,gscan);
image(I,'Tag','ToDelete');
axis image; 
g2scan = gscan + firstscan - 1;  
framename = deblank(F(g2scan,:));
t=title(framename,'FontSize',14,'Interpreter','none','Tag','ToDelete');
axis off;

%======================================================================== 
otherwise
%========================================================================
warning('Unknown action string')
%=======================================================================
end;
   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
