function [SAMimage, coords, stepsize] = ctf_read_svl(filename)
% ctf_read_svl - reads SAM image from a .svl file (generated by SAMsrc)
%
% [SAMimage, coords, stepsize] = ctf_read_svl(filename)
% 
% Returns SAMimage as voxels x 1 array of intensities
% coords contains corresponding voxel locations
%

% $Revision: 1.1 $ $Date: 2009-01-30 03:49:27 $

% Copyright (C) 2005, Sarang Dalal and Darren Weber
% 
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License
% as published by the Free Software Foundation; either version 2
% of the License, or (at your option) any later version.
% 
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
% 
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307,
% USA.

% History:  10/2005, Darren.Weber_at_radiology.ucsf.edu
%                    Sarang Dalal developed nut_read_svl.m
%                    DLW adapted nut_read_svl.m into ctf functions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ver = '$Revision: 1.1 $ $Date: 2009-01-30 03:49:27 $';
fprintf('CTF_READ_SVL [v %s]\n',ver(11:15));

if(nargin < 1)
    help ctf_read_svl;
    return;
end

fid=fopen(filename,'r','b');
identity = transpose(fread(fid,8,'*char'));
if(~strcmp(identity,'SAMIMAGE'))
    error('This doesn''t look like a SAM IMAGE file.');
end
vers = fread(fid,1,'int32');
setname = fread(fid,256,'*char');

% int32 vs. double for fiducial locations:
% The CTF svl format assumes that the fiducials will be 
% integer  values when defined in terms of MRI voxel 
% indices, but will be floating point values when 
% defined in MEG head coordinates.
numchans = fread(fid,1,'int32');
numweights = fread(fid,1,'int32');
if(numweights ~= 0)
    warning('... numweights ~= 0');
end

padbytes1 = fread(fid,1,'int32');

xstart = fread(fid,1,'double');
xend = fread(fid,1,'double');
ystart = fread(fid,1,'double');
yend = fread(fid,1,'double');
zstart = fread(fid,1,'double');
zend = fread(fid,1,'double');
stepsize = fread(fid,1,'double');

[meshz,meshy,meshx]=ndgrid(zstart:stepsize:zend,ystart:stepsize:yend,xstart:stepsize:xend);
coords = [meshx(:) meshy(:) meshz(:)];

hpfreq = fread(fid,1,'double');
lpfreq = fread(fid,1,'double');
bwfreq = fread(fid,1,'double');
meannoise = fread(fid,1,'double');

MRIname = transpose(fread(fid,256,'*char'));
nasion = fread(fid,3,'int32');
rightPA = fread(fid,3,'int32');
leftPA = fread(fid,3,'int32');

SAMtype = fread(fid,1,'int32');
SAMunit = fread(fid,1,'int32');

padbytes2 = fread(fid,1,'int32');

% DLW: why are these doubles and those above are int32?
nasion_meg = fread(fid,3,'double');
rightPA_meg = fread(fid,3,'double');
leftPA_meg = fread(fid,3,'double');

SAMunitname = fread(fid,32,'*char');

SAMimage = fread(fid,inf,'double');

fclose(fid);

return
