function values=extract(voxels,input_file_base,extension) 

%EXTRACT
%
% Extracts values of files at specified voxels.
%
% VALUES = EXTRACT ( VOXELS, INPUT_FILE_BASE, [EXTENSION] )
%
% VOXELS is a 3 column matrix of voxel locations in 'register' format
%
% INOUT_FILE_BASE is a matrix of file bases
%
% EXTENSION is a file extension to be added to INPUT_FILE_BASE.
% Default is [].
%
% VALUES is a #VOXELS x #FILEs matrix of extracted values.

if nargin<3
   extension=[]
end

numvoxels=size(voxels,1)
numfiles=size(input_file_base,1)
values=zeros(numvoxels,numfiles);
for file=1:numfiles
   in_handle=openimage([deblank(input_file_base(file,:)) extension]);
   if file==1
      numys=getimageinfo(in_handle, 'ImageHeight')
      numxs=getimageinfo(in_handle, 'ImageWidth')
      steps=getimageinfo(in_handle, 'Steps')
      i=voxels(:,1);
      if steps(1)<0
         i=numxs-1-i;
      end
      j=voxels(:,2);
      if steps(2)<0
         j=numys-1-j;
      end
      k=voxels(:,3);
   end
   for m=1:numvoxels
      X=reshape(getimages(in_handle,k(m)+1),numxs,numys);
      values(m,file)=X(i(m)+1,j(m)+1);
   end
end