function S = lyngby_getslice(V, plane, slice) 

% lyngby_getslice	- Extracts a slice from a volume
%
%	function S = lyngby_getslice(V, plane, slice) 
%
%	Input:	V	The volume in a vector
%		plane	's': sagittal, 'c': coronal, 't': transversal
%		slice	Index (index from one),	or 'mean' or 'max'
%
%	Output:	S	The slice as a matrix
%
%       Uses NUM_VOXELS, VOXEL_MASK and ROI_VOXELS to obtain the slice
%       from the vector.
%
%       See also: lyngby_global

% cvs : $Id: lyngby_getslice.m,v 1.7 1999/07/09 09:44:14 fnielsen Exp $ 
%       $Revision: 1.7 $

    lyngby_global;

    if length(V) == size(VOXEL_MASK, 2)
      % disp('irregular masked volume')
      V = V * VOXEL_MASK';
    end
    
    if (length(V) == prod(NUM_VOXELS))
      % No ROI_VOXELS
      dim = NUM_VOXELS;
      
      
      if plane == 's'
	VM = reshape(permute(reshape(V, dim), [2 3 1]), ...
	    prod(dim(2:3)), dim(1))'; 
	firstCoord = NUM_VOXELS(2);
	secondCoord = NUM_VOXELS(3); 
      elseif plane == 'c'
	VM = reshape(permute(reshape(V, dim), [1 3 2]), ...
	    prod(dim([1 3])), dim(2))'; 
	firstCoord = NUM_VOXELS(1);
	secondCoord = NUM_VOXELS(3); 
      elseif plane == 't'
	VM = reshape(V, prod(NUM_VOXELS(1:2)), NUM_VOXELS(3))';
	firstCoord = NUM_VOXELS(1);
	secondCoord = NUM_VOXELS(2); 
      else
	error('Wrong plane specification');
      end
      if isstr(slice)
	if size(VM,1) == 1
	  S = reshape(VM, firstCoord, secondCoord);
	elseif strcmp(slice, 'mean')
	  S = reshape(mean(VM), firstCoord, secondCoord);
	elseif strcmp(slice, 'max')
	  S = reshape(max(VM), firstCoord, secondCoord);
	else
	  error(sprintf('Slice ''%s'' specification not recognized')); 
	end
      else
	S = reshape(VM(slice, :),  firstCoord, secondCoord);
      end
    elseif ~isempty(ROI_VOXELS)
      % ROI_VOXELS mask
      
      dim = (ROI_VOXELS(:,2) - ROI_VOXELS(:,1) + 1)';
      if length(V) == prod(dim)
	%disp('ROIed volume')
	if plane == 's'
	  VM = reshape(permute(reshape(V, dim), [2 3 1]), ...
	      prod(dim(2:3)), dim(1))'; 
	  firstCoord = dim(2);
	  secondCoord = dim(3); 
	  slices = ROI_VOXELS(1,1):ROI_VOXELS(1,2);
	elseif plane == 'c'
	  VM = reshape(permute(reshape(V, dim), [1 3 2]), ...
	      prod(dim([1 3])), dim(2))'; 
	  firstCoord = dim(1);
	  secondCoord = dim(3); 
	  slices = ROI_VOXELS(2,1):ROI_VOXELS(2,2);
	elseif plane == 't'
	  VM = reshape(V, prod(dim(1:2)), dim(3))';
	  firstCoord = dim(1);
	  secondCoord = dim(2); 
	  slices = ROI_VOXELS(3,1):ROI_VOXELS(3,2);
	else
	  error('Wrong plane specification');
	end
	if isstr(slice)
	  if size(VM,1)==1
	    S = reshape(VM, firstCoord, secondCoord);
	  elseif strcmp(slice, 'mean')
	    S = reshape(mean(VM), firstCoord, secondCoord);
	  elseif strcmp(slice, 'max')
	    S = reshape(max(VM), firstCoord, secondCoord);
	  else
	    error(sprintf('Slice ''%s'' specification not recognized')); 
	  end
	else
	  s = find(slices == slice);
	  if isempty(s)
	    error(sprintf(['Specified slice is not in the volume. ' ...
		  'Slice: %d, Volume: %d-%d'], slice, ROI_VOXELS(3,1), ...
		ROI_VOXELS(3,2)));
	  else
	    S = reshape(VM(s, :), firstCoord, secondCoord);
	  end
	end
      else
	error('ROI_VOXELS and V are not of the same size.')
      end
    else
      error(['NUM_VOXELS, ROI_VOXELS, or perhaps VOXEL_MASK do not '...
	    'fit with the size of volume V.'])
    end





