function O = lyngby_read_volumes(filenames, varargin)

% lyngby_read_volumes   - Reads volumes from files
%
%	function O = lyngby_read_volumes(filename)
%
%	Input:    filename  Filenames as a cellstring
%
%       Property: Output       [ {mat} ]
%                 FilenameMask Filename for a file used as mask 
%                 VolMask      'Vol' structure with voxel mask
%
%	Output:   O            Output, eg, a 'mat' structure
%
%       Reads a number of ANALYZE files. The volumes are returned in
%       the 'matrix' field of a 'mat' structure with each row
%       corresponding to a volume and each column corresponding to the
%       time-series of a voxel. 
%
%       The voxel mask that is setup with the 'FilenameMask' and
%       'VolMask' property discards the voxels that a zero or NaN. If
%       the 'VolMask' is empty then all voxels will be loaded.
% 
%       See also LYNGBY, LYNGBY_READ_ANALYZE
%
% $Id: lyngby_read_volumes.m,v 1.3 2004/03/08 21:27:16 fnielsen Exp $ 


    % Check input arguments
    if nargin < 1
      error('Too few input arguments')
    end


    % Default properties
    filenameMask = ''; 
    indexMask    = [];
    orientation  = 'lrpasi';
    output       = 'mat';
    Vmask        = [];

    
    % Read properties
    for n = 1:2:length(varargin)
      arg = lower(varargin{n});

      if strcmp(arg, 'filenamemask') 
      	arg = varargin{n+1};
	if isstr(arg) 
	  filenameMask = arg;
	else
	  error('Argument with ''FilenameMask'' should be a string');
	end

      elseif strcmp(arg, 'volmask') 
      	arg = varargin{n+1};
	if isstruct(arg) | isempty(arg)
	  Vmask = arg;
	else
	  error('Argument with ''VolMask'' should be a structure');
	end
	
      else
	error(sprintf('Invalid property: %s', arg));
      end
    end

    
    % Voxel mask
    if isempty(Vmask)
      if ~isempty(filenameMask)
	siz = lyngby_read_analyze(filenameMask)
	v = lyngby_read_analyze(filenameMask);
	Vmask.volume = reshape(v, siz);
      end
    end
    siz = lyngby_read_analyze(filenames{1}, 'output', 'size');
    if isempty(Vmask)
      indexMask = 1:prod(siz);
    else
      indexMask = intersect(find(Vmask.volume(:)), find(~isnan(Vmask.volume)));
    end
    if isempty(indexMask)
      warning('No voxels in mask');
    end

    if strcmp(output, 'mat')
      
      % Allocate matrix
      O.matrix = zeros(length(filenames), length(indexMask));

      V.origin = lyngby_read_analyze(filenames{1}, 'output', ...
	  'origin');
      V.voxelSize = lyngby_read_analyze(filenames{1}, 'output', 'voxeldim');
      
      if strcmp('lr', orientation(1:2))
	V.x = ((1:siz(1)) - V.origin(1)) * V.voxelSize(1); 
      else
	V.x = ((1:siz(1)) - (siz(1)+1-V.origin(1))) * V.voxelSize(1);
      end
      if strcmp('pa', orientation(3:4))
	V.y = ((1:siz(2)) - V.origin(2)) * V.voxelSize(2);
      else
	V.y = ((1:siz(2)) - (siz(2)+1-V.origin(2))) * V.voxelSize(2);
      end
      if strcmp('is', orientation(5:6))
	V.z = ((1:siz(3)) - origin(3)) * V.voxelSize(3);
      else
	V.z = ((1:siz(3)) - (siz(3)+1-V.origin(3))) * V.voxelSize(3);
      end

      
      % Iterate over filenames
      if ~isempty(indexMask)
	for n = 1:length(filenames)
	  
	  % Read file
	  v = lyngby_read_analyze(filenames{n});
	  % Apply mask
	  O.matrix(n,:) = v(indexMask);
	end
      else
	% No voxel in mask
	O.matrix = zeros(length(filenames),0);
      end
      O.rows = filenames;
      O.type = 'mat';
      
    else
      error('Internal error')
    end



























