function X = lyngby_getdata(ROI, voxelMask, timeMask);

% lyngby_getdata       - Returns the masked datamatrix
%
%       function X = lyngby_getdata(ROI, voxelMask, timeMask);
%
%       Input:  ROI        (Optional) Region of Interests, 3x2 matrix
%               voxelMask  (Optional) Sparse masking matrix
%               timeMask   (Optional) Sparse masking matrix
%	
%       Output: X  The datamatrix
%
%       This function loads the datamatrix. The voxel is in the
%       horizontale direction (as columns) and the time is in the
%       vertical direction (as rows).
%       To save memory a 'Region of Interest' can be specified, 
%       and further a voxel or/and a time mask can be specified. 
%       If the input arguments are missing TIME_MASK, VOXEL_MASK and
%       ROI_VOXELS are used.
%
%       Which volumes is loaded is determined by global variables and
%       through calls to consecutive calls lyngby_getvolume (for each
%       volume). 
%
%       See also LYNGBY, LYNGBY_GETVOLUME, LYNGBY_UI_LOADFILE.
%
% $Id: lyngby_getdata.m,v 1.11 2002/05/19 22:27:00 fnielsen Exp $



    lyngby_global
        
    if nargin<3,
      timeMask = TIME_MASK;
      if nargin<2
	voxelMask = VOXEL_MASK;
	if nargin<1
	  ROI = ROI_VOXELS;
	end
      end
    end
    
    if issparse(timeMask)
      [fulltime,masktime] = find(timeMask);
    else
      masktime = (FILENAME_STARTINDEX:NUM_SCANS+FILENAME_STARTINDEX-1)';
    end

    if issparse(voxelMask)
      [maskvoxel, fullvoxel] = find(voxelMask);
    else
      if [3 2] == size(ROI);
	maskvoxel = 1:prod(ROI(:,2)-ROI(:,1)+1);
      else
	maskvoxel = 1:prod(NUM_VOXELS);
      end
    end
    
    rX = length(masktime);
    cX = length(maskvoxel);
    
    % Allocate big Datamatrix	
    lyngby_log('Allocating datamatrix ...')
    X = zeros(rX, cX);
    lyngby_log('Datamatrix allocated!')
	
    n = 1;
    for t = masktime',
      lyngby_log(sprintf('Reading volume number %i',t));		
      fullVolume = lyngby_getvolume(t);

      if [3 2] == size(ROI);
	x = ROI(1,1):ROI(1,2);
	y = ROI(2,1):ROI(2,2);
	z = ROI(3,1):ROI(3,2);
	yy = NUM_VOXELS(1)*(y-1);
	zz = prod(NUM_VOXELS(1:2))*(z-1);
	xy = kron(ones(size(y)),x)+kron(yy,ones(size(x)));
	xyz = kron(ones(size(zz)),xy)+kron(zz,ones(size(xy)));
	
	roiVolume = fullVolume(xyz);
	X(n,:) = roiVolume * voxelMask;
	
      else
	X(n,:) = fullVolume * voxelMask;
      end
      n = n + 1;
    end

 












