mne_load_coil_def

   [CoilDef,Header] = mne_load_coil_def(fname);
   CoilDef          = mne_load_coil_def(fname);

   If file name is not specified, the standard coil definition file
   $MNE_ROOT/setup/mne/coil_def.dat or $MNE_ROOT/share/mne/coil_def.dat is read

   The content of the coil definition file is described in
   section 5.6 of the MNE manual

   This routine is modified from the original BrainStorm routine
   created by John C. Mosher

0001 function [CoilDef,Header] = mne_load_coil_def(fname);
0002 %
0003 %
0004 %   [CoilDef,Header] = mne_load_coil_def(fname);
0005 %   CoilDef          = mne_load_coil_def(fname);
0006 %
0007 %   If file name is not specified, the standard coil definition file
0008 %   $MNE_ROOT/setup/mne/coil_def.dat or $MNE_ROOT/share/mne/coil_def.dat is read
0009 %
0010 %   The content of the coil definition file is described in
0011 %   section 5.6 of the MNE manual
0012 %
0013 %   This routine is modified from the original BrainStorm routine
0014 %   created by John C. Mosher
0015 %
0016 
0017 %
0018 %   John C. Mosher
0019 %   Los Alamos National Laboratory
0020 %
0021 %   Author : Matti Hamalainen, MGH Martinos Center
0022 %   License : BSD 3-clause
0023 %
0024 %   Copyright (c) 2005 BrainStorm MMIV by the University of Southern California
0025 %   Principal Investigator:
0026 %   ** Professor Richard M. Leahy, USC Signal & Image Processing Institute
0027 %
0028 %
0029 %   Revision 1.5  2006/09/08 19:27:13  msh
0030 %   Added KIT coil type to mne_load_coil_def
0031 %   Allow reading of measurement info by specifying just a file name.
0032 %
0033 %   Revision 1.4  2006/04/25 12:29:10  msh
0034 %   Added Magnes and CTF reference coil drawings.
0035 %
0036 %   Revision 1.3  2006/04/23 15:29:40  msh
0037 %   Added MGH to the copyright
0038 %
0039 %   Revision 1.2  2006/04/17 15:01:34  msh
0040 %   More small improvements.
0041 %
0042 %   Revision 1.1  2006/04/17 11:52:15  msh
0043 %   Added coil definition stuff
0044 %
0045 %
0046 
0047 me='MNE:mne_load_coil_def';
0048 
0049 %% Setup the default inputs
0050 
0051 if nargin == 0
0052     fname = mne_file_name('setup/mne/coil_def.dat');
0053     if ~exist(fname,'file')
0054         fname = mne_file_name('share/mne/coil_def.dat');
0055         if ~exist(fname,'file')
0056             error(me,'The standard coil definition file was not found');
0057         end
0058     end
0059 end
0060 
0061 % Read in the coil_def information
0062 
0063 %% Open, Read in entire file, close
0064 
0065 fid = fopen(fname,'rt');
0066 
0067 if fid < 0
0068     error(me,'Could not open coil definition file %s',fname);
0069 end
0070 
0071 istr = 1; % string indexer
0072 str_array = cell(1000,1); % preallocate
0073 str_array{istr} = fgetl(fid); % get first string
0074 
0075 while ischar(str_array{istr}),
0076     istr = istr + 1;
0077     str_array{istr} = fgetl(fid); % get next string
0078 end
0079 fclose(fid);
0080 
0081 str_array = str_array(1:(istr-1)); % trim allocation
0082 
0083 %% Read the Header, find the structure
0084 
0085 % Gather the header lines
0086 HeaderLines = strmatch('#',str_array); % lines that begin with a comment
0087 
0088 % where are the structure lines
0089 StructureLines = strmatch('# struct',str_array); % subset of header lines
0090 % should be two
0091 if length(StructureLines) ~= 2,
0092     error(me,'%s anticipates two lines of structures',mfilename)
0093 end
0094 
0095 % with each structure line is a format line
0096 FormatLines = strmatch('# format',str_array); % subset of header lines
0097 % assume there are also two
0098 
0099 FieldNames = cell(1,2);
0100 Format = cell(1,2);
0101 % first structure is the coil information
0102 % won't actually use the second structure, just its format
0103 for i = 1:2,
0104     FieldNames{i}  = strread(str_array{StructureLines(i)},'%s');
0105     FieldNames{i}(1:2) = []; % strip the comment symbol and struct keyword
0106     
0107     [ignore,Format{i}] = strtok(str_array{FormatLines(i)},'''');
0108     Format{i} = strrep(Format{i},'''',''); % strip the single quotes
0109 end
0110 
0111 %% Allocate the arrays for loading
0112 
0113 % interleave every fieldname with a null value
0114 struct_arg = [FieldNames{1} cell(length(FieldNames{1}),1)]';
0115 % each column an argument pair
0116 
0117 % Preallocate a structure
0118 [CoilDef(1:100)] = deal(struct(struct_arg{:}));
0119 
0120 
0121 % Convert the rest of the string array to a structure
0122 iCoil = 0; % counter
0123 iLine = HeaderLines(end); % next line after header
0124 
0125 while iLine < length(str_array), % number of lines in file
0126     iCoil = iCoil + 1; % next coil definition
0127     iLine = iLine + 1; % next line
0128     
0129     % first read the integer information on the coil
0130     % begin by breaking the line into two parts, numeric and description
0131     [numeric_items, description] = strtok(str_array{iLine},'"');
0132     temp = sscanf(numeric_items,Format{1}); % extra %s doesn't matter
0133     % assign temp in the order of the fields
0134     for i = 1:(length(FieldNames{1})-1),
0135         CoilDef(iCoil).(FieldNames{1}{i}) = temp(i);
0136     end
0137     % then assign the description
0138     % let's strip the quotes first
0139     description = strrep(description,'"','');
0140     CoilDef(iCoil).(FieldNames{1}{end}) = description;
0141     
0142     % now read the coil definition
0143     CoilDef(iCoil).coildefs = zeros(CoilDef(iCoil).num_points,7);
0144     for i = 1:CoilDef(iCoil).num_points,
0145         iLine = iLine + 1;
0146         CoilDef(iCoil).coildefs(i,:) = sscanf(str_array{iLine},...
0147             Format{2})';
0148     end
0149     
0150     % now draw it
0151     % local subfunction below
0152     CoilDef(iCoil).FV = draw_sensor(CoilDef(iCoil));
0153     
0154 end
0155 
0156 CoilDef = CoilDef(1:iCoil); % trim allocation
0157 
0158 Header = str_array(HeaderLines);
0159 
0160 function FV = draw_sensor(CoilDef);
0161 
0162 % create a patch based on the sensor type
0163 
0164 % The definitions as of 14 October 2005:
0165 % for i = 1:3:length(CoilDef),fprintf('%d %s\n',CoilDef(i).id,CoilDef(i).description);end
0166 % 2 Neuromag-122 planar gradiometer size = 27.89  mm base = 16.20  mm
0167 % 2000 Point magnetometer
0168 % 3012 Vectorview planar gradiometer T1 size = 26.39  mm base = 16.80  mm
0169 % 3013 Vectorview planar gradiometer T2 size = 26.39  mm base = 16.80  mm
0170 % 3022 Vectorview magnetometer T1 size = 25.80  mm
0171 % 3023 Vectorview magnetometer T2 size = 25.80  mm
0172 % 3024 Vectorview magnetometer T3 size = 21.00  mm
0173 % 4001 Magnes WH2500 magnetometer size = 11.50  mm
0174 % 4002 Magnes WH3600 gradiometer size = 18.00  mm base = 50.00  mm
0175 % 5001 CTF axial gradiometer size = 18.00  mm base = 50.00  mm
0176 % 7001 BabySQUID I axial gradiometer size = 6.0 mm base = 50 mm
0177 
0178 FV = struct('faces',[],'vertices',[]); % standard convention
0179 
0180 % recall that vertices are 1 per ROW, not column, of matrix
0181 
0182 switch CoilDef.id
0183     case {2,3012,3013,3011}
0184         % square figure eight
0185         % wound by right hand rule such that +x side is "up" (+z)
0186         LongSide = CoilDef.size*1000; % length of long side in mm
0187         Offset = 2.5; % mm offset of the center portion of planar grad coil
0188         FV.vertices = [0 0 0; Offset 0 0; ...
0189             Offset -LongSide/2 0; LongSide/2 -LongSide/2 0; ...
0190             LongSide/2 LongSide/2 0; ...
0191             Offset LongSide/2 0; Offset 0 0; ...
0192             0 0 0; -Offset 0 0; -Offset -LongSide/2 0; ...
0193             -LongSide/2 -LongSide/2 0; ...
0194             -LongSide/2 LongSide/2 0; ...
0195             -Offset LongSide/2 0; -Offset 0 0]/1000;
0196         FV.faces = [1:length(FV.vertices)];
0197     case 2000
0198         % point source
0199         LongSide = 2; % mm, tiny square
0200         FV.vertices = [-1 1 0;1 1 0;1 -1 0; -1 -1 0]*LongSide/1000/2;
0201         FV.faces = [1:length(FV.vertices)];
0202     case {3022, 3023, 3024}
0203         % square magnetometer
0204         LongSide = CoilDef.size*1000; % mm, length of one side
0205         FV.vertices = [-1 1 0;1 1 0;1 -1 0; -1 -1 0]*LongSide/1000/2;
0206         FV.faces = [1:length(FV.vertices)];
0207     case {4001,4003,5002}
0208         % round magnetometer
0209         Radius = CoilDef.size*1000/2; % mm, radius of coil
0210         Len_cir = 15; % number of points for circle
0211         circle = cos(2*pi*[0:(Len_cir-1)]/Len_cir) + ...
0212             sqrt(-1)*sin(2*pi*[0:(Len_cir-1)]/Len_cir); % complex circle unit
0213         FV.vertices = ...
0214             [real(circle)' imag(circle)' zeros(Len_cir,1)]*Radius/1000;
0215         FV.faces = [1:length(FV.vertices)];
0216     case {4002, 5001, 5003, 4004, 6001, 7001}
0217         % round coil 1st order gradiometer
0218         Radius = CoilDef.size*1000/2; % mm radius
0219         Baseline = CoilDef.baseline*1000; % axial separation
0220         Len_cir = 15; % number of points for circle
0221         % This time, go all the way around circle to close it fully
0222         circle = cos(2*pi*[0:Len_cir]/Len_cir) + ...
0223             sqrt(-1)*sin(2*pi*[0:Len_cir]/Len_cir); % complex circle unit
0224         circle = circle*Radius; % scaled
0225         FV.vertices = ...
0226             [[real(circle)' imag(circle)' zeros(Len_cir+1,1)];... % first coil
0227             [real(circle)' -imag(circle)' zeros(Len_cir+1,1)+Baseline]]/1000; % 2nd coil
0228         FV.faces = [1:length(FV.vertices)];
0229     case {5004,4005}
0230         % round coil 1st order off-diagonal gradiometer
0231         Radius = CoilDef.size*1000/2; % mm radius
0232         Baseline = CoilDef.baseline*1000; % axial separation
0233         Len_cir = 15; % number of points for circle
0234         % This time, go all the way around circle to close it fully
0235         circle = cos(2*pi*[0:Len_cir]/Len_cir) + ...
0236             sqrt(-1)*sin(2*pi*[0:Len_cir]/Len_cir); % complex circle unit
0237         circle = circle*Radius; % scaled
0238         FV.vertices = ...
0239             [[real(circle)'+Baseline/2.0 imag(circle)' zeros(Len_cir+1,1)];... % first coil
0240             [real(circle)'-Baseline/2.0 -imag(circle)' zeros(Len_cir+1,1)]]/1000; % 2nd coil
0241         FV.faces = [1:length(FV.vertices)];
0242     otherwise
0243         FV = [];
0244 end
0245 
0246