read_ctf_sens.m

function [magn] = read_ctf_sens(filename);

% READ_CTF_SENS reads MEG sensor information from CTF configuration file
%
% magn = read_ctf_sens(filename)
% 
% where the returned structure meg has the fields
%   magn.pnt	position first coil
%   magn.ori	orientation first coil
%   magn.pnt2	position second coil
%   magn.ori2	orientation second coil

% Copyright (C) 2002, Robert Oostenveld
% 
fid = fopen(filename, 'r');
if fid==-1
  error(sprintf('could not open file %s', filename));
end

% skip the first line
line = fgetl(fid);

% read the channel parameters
chan = 0;
while (1)
  line = fgetl(fid);
  if ~isempty(line) & line==-1
    % reached end of file
    break
  end

  element = tokenize(line);
  if ~isempty(element)
    chan = chan+1;

    % channel label or name
    name{chan} = strtok(element{1}, '-');

    % type of channel (EEG/MEG etc.), convert to uppercase
    type{chan} = upper(element{7});

    % number of coils (magnetometer=1/gradiometer=2)
    ncoils(chan) = str2num(element{2});

    % polarity of channel
    polarity(chan) = str2num(element{26});

    % position of first coil
    coil1.pnt(chan,:) = [str2num(element{17}) str2num(element{18}) str2num(element{19})];
    % direction of first coil
    coil1.dir(chan,:) = [str2num(element{20}) str2num(element{21})];

    % length of basline
    baseline.length(chan) = str2num(element{16});
    % direction of baseline
    baseline.dir(chan,:) = [str2num(element{22}) str2num(element{23})];

    % direction of second coil
    coil2.dir(chan,:) = [str2num(element{24}) str2num(element{25})];
  end
end

% find the sensor gradiometer channels
indx = find(strcmp(type, 'GRAD1-SENS'));
nchan = length(indx);
label = name(indx);
pnt1 = coil1.pnt(indx,:) * 10;	% convert from cm to mm

% shift offset to origin to make triangulation
offset = [mean(pnt1(:,1)) mean(pnt1(:,2)) min(pnt1(:,3))];
prj = elproj(pnt1 - repmat(offset, nchan, 1));
tri = delaunay(prj(:,1), prj(:,2));

% compute the directional vectors of each gradiometer: bottom coil
a = coil1.dir(indx,1);
b = coil1.dir(indx,2);
[r1, r2, r3] = sph2cart(a*2*pi/360, pi/2 - b*2*pi/360, ones(nchan,1));
r1 = -r1 .* polarity(indx)';
r2 = -r2 .* polarity(indx)';
r3 = -r3 .* polarity(indx)';
ori1 = [r1 r2 r3];

% compute the directional vectors of each gradiometer: top coil
a = coil2.dir(indx,1);
b = coil2.dir(indx,2);
[r1, r2, r3] = sph2cart(a*2*pi/360, pi/2 - b*2*pi/360, ones(nchan,1));
r1 = -r1 .* polarity(indx)';
r2 = -r2 .* polarity(indx)';
r3 = -r3 .* polarity(indx)';
ori2 = [r1 r2 r3];

% compute the directional vectors of the baseline and the position of the 2nd coil
a = baseline.dir(indx,1);
b = baseline.dir(indx,2);
[r1, r2, r3] = sph2cart(a*2*pi/360, pi/2 - b*2*pi/360, ones(nchan,1));
oriB = [r1 r2 r3];
pnt2 = pnt1 + oriB .* repmat(baseline.length(indx)', 1, 3);

magn.pnt	= pnt1;
magn.ori	= ori1;
magn.pnt2	= pnt2;
magn.ori2	= ori2;
magn.tri	= tri;
magn.label	= label;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [element] = tokenize(line);
count = 0;
while ~isempty(line)
  [tok, line] = strtok(line);
  if ~isempty(tok)
    count = count+1;
    element{count} = tok;
  else
    break;
  end
end