mne_add_coil_defs

 [res] = mne_add_coil_defs(chs,accuracy,coil_def_templates)

 Add transformed coil definitions to the channel info

 chs        - original channel definitions
 accuracy   - desired accuracy (0, 1, or 2, defaults to 1)
 templates  - coil definition templates
              (defaults to $MNE_ROOT/setup/mne/coil_def.dat or $MNE_ROOT/share/mne/coil_def.dat)

0001 function [res] = mne_add_coil_defs(chs,accuracy,templates)
0002 %
0003 % [res] = mne_add_coil_defs(chs,accuracy,coil_def_templates)
0004 %
0005 % Add transformed coil definitions to the channel info
0006 %
0007 % chs        - original channel definitions
0008 % accuracy   - desired accuracy (0, 1, or 2, defaults to 1)
0009 % templates  - coil definition templates
0010 %              (defaults to $MNE_ROOT/setup/mne/coil_def.dat or $MNE_ROOT/share/mne/coil_def.dat)
0011 %
0012 
0013 %
0014 %   Author : Matti Hamalainen, MGH Martinos Center
0015 %   License : BSD 3-clause
0016 %
0017 %
0018 %   Revision 1.8  2006/05/03 18:53:05  msh
0019 %   Approaching Matlab 6.5 backward compatibility
0020 %
0021 %   Revision 1.7  2006/04/26 00:14:24  msh
0022 %   Return empty structure from fiff_pick_channels evoked and fiff_pick_types_evoked if there is no match.
0023 %   Accuracy checking was incorrect in mne_add_coil_defs
0024 %
0025 %   Revision 1.6  2006/04/24 00:07:24  msh
0026 %   Fixed error in the computation of the surface-based coordinate system.
0027 %
0028 %   Revision 1.5  2006/04/23 15:29:40  msh
0029 %   Added MGH to the copyright
0030 %
0031 %   Revision 1.4  2006/04/21 17:31:07  msh
0032 %   Added the examples.
0033 %   Modified the formatting of some informative output.
0034 %
0035 %   Revision 1.3  2006/04/18 20:44:46  msh
0036 %   Added reading of forward solution.
0037 %   Use length instead of size when appropriate
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 me='MNE:mne_add_coil_defs';
0047 
0048 global FIFF;
0049 if isempty(FIFF)
0050     FIFF = fiff_define_constants();
0051 end
0052 
0053 if nargin == 2 || nargin == 1 || isempty(templates)
0054     %
0055     %   We need the templates
0056     %
0057     try
0058         templates = mne_load_coil_def;
0059     catch
0060         error(me,'%s',mne_omit_first_line(lasterr));
0061     end
0062     if nargin == 1
0063         accuracy = 1;
0064     end
0065 elseif nargin ~= 3
0066     error(me,'Incorrect number of arguments');
0067 end
0068 
0069 fprintf(1,'\t%d coil definition templates available\n',length(templates));
0070 
0071 if accuracy ~= 0 && accuracy ~= 1 && accuracy ~= 2
0072     error(me,'Accuracy should attain one of the values 0, 1, 2');
0073 end
0074 
0075 res    = chs;
0076 nchan  = length(chs);
0077 ntemp  = length(templates);
0078 nmeg   = 0;
0079 neeg   = 0;
0080 eeg_cm = zeros(1,3);
0081 
0082 first = true;
0083 for k = 2:nchan
0084     if chs(k).kind == FIFF.FIFFV_EEG_CH || chs(k).kind == FIFF.FIFFV_MEG_CH || ...
0085             chs(k).kind == FIFF.FIFFV_REF_MEG_CH
0086         if first
0087             coord_frame = chs(k).coord_frame;
0088             first = false;
0089         elseif chs(k).coord_frame ~= coord_frame
0090             error(me,'All coils and electrodes are not in the same coordinate system (%s)',chs(k).ch_name);
0091         end
0092     end
0093 end
0094 
0095 lower_half = hex2dec('FFFF');
0096 for k = 1:nchan
0097     res(k).coil_def = [];
0098     if chs(k).kind == FIFF.FIFFV_EEG_CH && ~isempty(chs(k).eeg_loc)
0099         %
0100         %   Set up a "coil definition" for an EEG electrode
0101         %
0102         temp             = templates(1);             %   Just pick one
0103         temp.class       = 1000;                     %   EEG electrode class
0104         temp.id          = chs(k).coil_type;         %   Copy from the channel info
0105         temp.accuracy    = 1;                        %   Normal accuracy no matter what
0106         temp.num_points  = size(chs(k).eeg_loc,2);   %   Number of points is normally two
0107         temp.size        = 10e-3;                    %   Rather arbitrary
0108         temp.baseline    = 0;                        %   Baseline does not apply
0109         temp.description = 'EEG electrode';          %   Useful description
0110         temp.coildefs    = zeros(temp.num_points,7); %   Location data
0111         temp.coildefs(1,1)   = 1;                    %   The electrode location
0112         temp.coildefs(1,2:4) = chs(k).eeg_loc(:,1)';
0113         nn  = chs(k).eeg_loc(:,1);                   %   The normal definition is probably immaterial
0114         len = nn'*nn;
0115         if len > 0
0116             temp.coildefs(1,5:7) = nn';
0117         end
0118         if temp.num_points > 1                       %   Add the reference electrode
0119             temp.coildefs(2,1)   = -1;
0120             temp.coildefs(2,2:4) = chs(k).eeg_loc(:,2);
0121             nn  = chs(k).eeg_loc(:,2);
0122             len = nn'*nn;
0123             if len > 0
0124                 temp.coildefs(2,5:7) = nn';
0125             end
0126         end
0127         temp.FV          = [];                            %   Empty for now
0128         eeg_cm           = eeg_cm + temp.coildefs(1,2:4); %   Accumulate center of mass
0129         res(k).coil_def  = temp;                          %   There it is
0130         neeg = neeg + 1;
0131     elseif chs(k).kind == FIFF.FIFFV_MEG_CH || chs(k).kind == FIFF.FIFFV_REF_MEG_CH
0132         %
0133         %   Set up an MEG coil definition
0134         %
0135         coil_type = bitand(double(chs(k).coil_type),lower_half);
0136         temp = [];
0137         for t = 1:ntemp
0138             if templates(t).id == coil_type && templates(t).accuracy == accuracy
0139                 temp = templates(t);
0140                 break;
0141             end
0142         end
0143         %
0144         %  Did we find the template?
0145         %
0146         if isempty(temp)
0147             error(me,'Could not find an MEG coil template (coil type = %d accuracy = %d) for channel %s', ...
0148                 coil_type,accuracy,chs(k).ch_name);
0149         end
0150         %
0151         %  Transform the template using the coil transformation
0152         %
0153         np = size(temp.coildefs,1);
0154         %
0155         %  Transform the integration points and the normals
0156         %
0157         trans = chs(k).coil_trans(1:3,:);
0158         temp.coildefs(:,2:4) = (trans*[ temp.coildefs(:,2:4)' ; ones(1,np) ])';
0159         temp.coildefs(:,5:7) = (trans*[ temp.coildefs(:,5:7)' ; zeros(1,np) ])';
0160         %
0161         %  Transform the vertices
0162         %
0163         if ~isempty(temp.FV)
0164             nvert = size(temp.FV.vertices,1);
0165             temp.FV.vertices = (trans*[temp.FV.vertices' ; ones(1,nvert)])';
0166         end
0167         %
0168         %   Attach the coil definition to the result
0169         %
0170         res(k).coil_def = temp;
0171         nmeg = nmeg + 1;
0172     end
0173 end
0174 %
0175 %     Create the EEG coils drawings now that we have the center of mass
0176 %
0177 if neeg > 0
0178     eeg_cm = eeg_cm/neeg;
0179     for k = 1:nchan
0180         if res(k).kind == FIFF.FIFFV_EEG_CH && ~isempty(res(k).eeg_loc)
0181             res(k).coil_def.FV = draw_electrode(res(k).coil_def,eeg_cm);
0182         end
0183     end
0184 end
0185 
0186 fprintf(1,'\t%d MEG coil definitions and %d EEG electrodes set up\n',nmeg,neeg);
0187 
0188 return;
0189 
0190     function FV = draw_electrode(def,cm)
0191 
0192         % create a patch for an EEG electrode
0193 
0194         if def.class == 1000
0195             % round disk
0196             Radius  = def.size/2;           % radius
0197             Len_cir = 15;                   % number of points for circle
0198             r0      = def.coildefs(1,2:4);
0199             %
0200             %   Make a template circle
0201             %
0202             rr(1,:) = cos(2*pi*[0:Len_cir]/Len_cir);
0203             rr(2,:) = sin(2*pi*[0:Len_cir]/Len_cir);
0204             rr(3,:) = 0.0;
0205             rr(4,:) = ones(1,Len_cir+1);
0206             %
0207             %    Scaling
0208             %
0209             T1 = [ [ Radius 0 0 0 ] ; [ 0 Radius 0 0 ] ; [ 0 0 1 0 ...
0210                 ] ; [  0 0 0 1 ] ];
0211             %
0212             %   Rotation (this is only approximately correct...)
0213             %
0214             %   The normal is parallel to the vector from CM to the electrode location
0215             %
0216             nn     = (def.coildefs(1,2:4) - cm)';
0217             sizenn = norm(nn);
0218             if sizenn > 0
0219                 nn = nn/sizenn;
0220                 [ u, s, v ]  = svd(eye(3,3) - nn*nn');
0221                 if nn'*u(:,3) < 0
0222                     u = - u;
0223                 end
0224                 T2 = [ u [ 0 0 0 ]' ; [ 0 0 0 1 ]];
0225             else
0226                 T2 = eye(4,4);
0227             end
0228             %
0229             %   Translation
0230             %
0231             T3 = [ [ 1 0 0 r0(1) ] ; [ 0 1 0 r0(2) ] ; [ 0 0 1 r0(3) ...
0232                 ] ; [  0 0 0 1 ] ];
0233             %
0234             %    Transfrom the data and create the structure
0235             %
0236             rr = T3*T2*T1*rr;
0237             FV.vertices = [[rr(1,:)' rr(2,:)' rr(3,:)']];
0238             FV.faces    = [1:length(FV.vertices)];
0239         end
0240         return;
0241     end
0242 
0243 end