mne_read_forward_solution

 [fwd] = mne_read_forward_solution(fname,force_fixed,surf_ori,include,exclude)

 A forward solution from a fif file

 fname        - The name of the file
 force_fixed  - Force fixed source orientation mode? (optional)
 surf_ori     - Use surface based source coordinate system? (optional)
 include      - Include these channels (optional)
 exclude      - Exclude these channels (optional)

0001 function [fwd] = mne_read_forward_solution(fname,force_fixed,surf_ori,include,exclude)
0002 %
0003 % [fwd] = mne_read_forward_solution(fname,force_fixed,surf_ori,include,exclude)
0004 %
0005 % A forward solution from a fif file
0006 %
0007 % fname        - The name of the file
0008 % force_fixed  - Force fixed source orientation mode? (optional)
0009 % surf_ori     - Use surface based source coordinate system? (optional)
0010 % include      - Include these channels (optional)
0011 % exclude      - Exclude these channels (optional)
0012 %
0013 
0014 %
0015 %   Author : Matti Hamalainen, MGH Martinos Center
0016 %   License : BSD 3-clause
0017 %
0018 %   Revision 1.11  2008/06/07 21:22:10  msh
0019 %   Added FIFF_FORWARD_SOLUTION_GRAD constant
0020 %
0021 %   Revision 1.10  2008/06/07 21:19:00  msh
0022 %   Added reading of the field gradients.
0023 %
0024 %   Revision 1.9  2006/11/30 05:43:29  msh
0025 %   Fixed help text in fiff_dir_tree_find
0026 %   Fixed check for the existence of parent MRI block in mne_read_forward_solution
0027 %
0028 %   Revision 1.8  2006/05/05 03:50:40  msh
0029 %   Added routines to compute L2-norm inverse solutions.
0030 %   Added mne_write_inverse_sol_stc to write them in stc files
0031 %   Several bug fixes in other files
0032 %
0033 %   Revision 1.7  2006/05/03 19:09:03  msh
0034 %   Fixed even more compatibility issues.
0035 %
0036 %   Revision 1.6  2006/04/24 00:07:24  msh
0037 %   Fixed error in the computation of the surface-based coordinate system.
0038 %
0039 %   Revision 1.5  2006/04/23 15:29:41  msh
0040 %   Added MGH to the copyright
0041 %
0042 %   Revision 1.4  2006/04/21 21:33:07  msh
0043 %   Fixed error in combination of the forward solution.
0044 %   Fixed help text in mne_read_inverse_operator
0045 %
0046 %   Revision 1.3  2006/04/20 21:49:38  msh
0047 %   Added mne_read_inverse_operator
0048 %   Changed some of the routines accordingly for more flexibility.
0049 %
0050 %   Revision 1.2  2006/04/18 23:21:22  msh
0051 %   Added mne_transform_source_space_to and mne_transpose_named_matrix
0052 %
0053 %   Revision 1.1  2006/04/18 20:44:46  msh
0054 %   Added reading of forward solution.
0055 %   Use length instead of size when appropriate
0056 %
0057 %
0058 
0059 me='MNE:mne_read_forward_solution';
0060 
0061 global FIFF;
0062 if isempty(FIFF)
0063     FIFF = fiff_define_constants();
0064 end
0065 
0066 if nargin == 1
0067     include = [];
0068     exclude = [];
0069     force_fixed = false;
0070     surf_ori = false;
0071 elseif nargin == 2
0072     include = [];
0073     exclude = [];
0074     surf_ori = false;
0075 elseif nargin == 3
0076     include = [];
0077     exclude = [];
0078 elseif nargin == 4
0079     exclude = [];
0080 elseif nargin ~= 5
0081     error(me,'Incorrect number of arguments');
0082 end
0083 %
0084 %   Open the file, create directory
0085 %
0086 fprintf(1,'Reading forward solution from %s...\n',fname);
0087 [ fid, tree ] = fiff_open(fname);
0088 %
0089 %   Find all forward solutions
0090 %
0091 fwds = fiff_dir_tree_find(tree,FIFF.FIFFB_MNE_FORWARD_SOLUTION);
0092 if isempty(fwds)
0093     fclose(fid);
0094     error(me,'No forward solutions in %s',fname);
0095 end
0096 %
0097 %   Parent MRI data
0098 %
0099 parent_mri = fiff_dir_tree_find(tree,FIFF.FIFFB_MNE_PARENT_MRI_FILE);
0100 if isempty(parent_mri)
0101     fclose(fid);
0102     error(me,'No parent MRI information in %s',fname);
0103 end
0104 %
0105 %   Parent MEG data
0106 %
0107 parent_meg = fiff_dir_tree_find(tree,FIFF.FIFFB_MNE_PARENT_MEAS_FILE);
0108 if isempty(parent_meg)
0109     fclose(fid);
0110     error(me,'No parent MEG information in %s',fname);
0111 end
0112 p = 0;
0113 for k = 1:parent_meg.nent
0114     kind = parent_meg.dir(k).kind;
0115     pos  = parent_meg.dir(k).pos;
0116     if kind == FIFF.FIFF_CH_INFO
0117         p = p+1;
0118         tag = fiff_read_tag(fid,pos);
0119         chs(p) = tag.data;
0120     end
0121 end
0122 
0123 try
0124     src = mne_read_source_spaces(fid,false,tree);
0125 catch
0126     fclose(fid);
0127     error(me,'Could not read the source spaces (%s)',mne_omit_first_line(lasterr));
0128 end
0129 for k = 1:length(src)
0130     src(k).id = mne_find_source_space_hemi(src(k));
0131 end
0132 fwd = [];
0133 %
0134 %   Locate and read the forward solutions
0135 %
0136 megnode = [];
0137 eegnode = [];
0138 for k = 1:length(fwds)
0139     tag = find_tag(fwds(k),FIFF.FIFF_MNE_INCLUDED_METHODS);
0140     if isempty(tag)
0141         fclose(fid);
0142         error(me,'Methods not listed for one of the forward solutions');
0143     end
0144     if tag.data == FIFF.FIFFV_MNE_MEG
0145         megnode = fwds(k);
0146     elseif tag.data == FIFF.FIFFV_MNE_EEG
0147         eegnode = fwds(k);
0148     end
0149 end
0150 
0151 megfwd = read_one(megnode);
0152 if ~isempty(megfwd)
0153     if megfwd.source_ori == FIFF.FIFFV_MNE_FIXED_ORI
0154         ori = 'fixed';
0155     else
0156         ori = 'free';
0157     end
0158     fprintf(1,'\tRead MEG forward solution (%d sources, %d channels, %s orientations)\n',...
0159         megfwd.nsource,megfwd.nchan,ori);
0160 end
0161 eegfwd = read_one(eegnode);
0162 if ~isempty(eegfwd)
0163     if eegfwd.source_ori == FIFF.FIFFV_MNE_FIXED_ORI
0164         ori = 'fixed';
0165     else
0166         ori = 'free';
0167     end
0168     fprintf(1,'\tRead EEG forward solution (%d sources, %d channels, %s orientations)\n',...
0169         eegfwd.nsource,eegfwd.nchan,ori);
0170 end
0171 %
0172 %   Merge the MEG and EEG solutions together
0173 %
0174 if ~isempty(megfwd) && ~isempty(eegfwd)
0175     if size(megfwd.sol.data,2) ~= size(eegfwd.sol.data,2) || ...
0176             megfwd.source_ori ~= eegfwd.source_ori || ...
0177             megfwd.nsource ~= eegfwd.nsource || ...
0178             megfwd.coord_frame ~= eegfwd.coord_frame
0179         fclose(fid);
0180         error(me,'The MEG and EEG forward solutions do not match');
0181     end
0182     fwd = megfwd;
0183     fwd.sol.data      = [ fwd.sol.data ; eegfwd.sol.data ];
0184     fwd.sol.nrow      = fwd.sol.nrow + eegfwd.sol.nrow;
0185     fwd.sol.row_names = [ fwd.sol.row_names eegfwd.sol.row_names ];
0186     if  ~isempty(fwd.sol_grad)
0187         fwd.sol_grad.data      = [ fwd.sol_grad.data ; eegfwd.sol_grad.data ];
0188         fwd.sol_grad.nrow      = fwd.sol_grad.nrow + eegfwd.sol_grad.nrow;
0189         fwd.sol_grad.row_names = [ fwd.sol_grad.row_names eegfwd.sol_grad.row_names ];
0190     end
0191     fwd.nchan         = fwd.nchan + eegfwd.nchan;
0192     fprintf(1,'\tMEG and EEG forward solutions combined\n');
0193 elseif ~isempty(megfwd)
0194     fwd = megfwd;
0195 else
0196     fwd = eegfwd;
0197 end
0198 clear('megfwd');
0199 clear('eegfwd');
0200 %
0201 %   Get the MRI <-> head coordinate transformation
0202 %
0203 tag = find_tag(parent_mri,FIFF.FIFF_COORD_TRANS);
0204 if isempty(tag)
0205     fclose(fid);
0206     error(me,'MRI/head coordinate transformation not found');
0207 else
0208     mri_head_t = tag.data;
0209     if mri_head_t.from ~= FIFF.FIFFV_COORD_MRI || mri_head_t.to ~= FIFF.FIFFV_COORD_HEAD
0210         mri_head_t = fiff_invert_transform(mri_head_t);
0211         if mri_head_t.from ~= FIFF.FIFFV_COORD_MRI || mri_head_t.to ~= FIFF.FIFFV_COORD_HEAD
0212             fclose(fid);
0213             error(me,'MRI/head coordinate transformation not found');
0214         end
0215     end
0216 end
0217 fclose(fid);
0218 %
0219 fwd.mri_head_t = mri_head_t;
0220 %
0221 %   Transform the source spaces to the correct coordinate frame
0222 %   if necessary
0223 %
0224 if fwd.coord_frame ~= FIFF.FIFFV_COORD_MRI && ...
0225         fwd.coord_frame ~= FIFF.FIFFV_COORD_HEAD
0226     error(me,'Only forward solutions computed in MRI or head coordinates are acceptable');
0227 end
0228 %
0229 nuse = 0;
0230 for k = 1:length(src)
0231     try
0232         src(k) = mne_transform_source_space_to(src(k),fwd.coord_frame,mri_head_t);
0233     catch
0234         error(me,'Could not transform source space (%s)',mne_omit_first_line(lasterr));
0235     end
0236     nuse = nuse + src(k).nuse;
0237 end
0238 if nuse ~= fwd.nsource
0239     error(me,'Source spaces do not match the forward solution.\n');
0240 end
0241 fprintf(1,'\tSource spaces transformed to the forward solution coordinate frame\n');
0242 fwd.src = src;
0243 %
0244 %   Handle the source locations and orientations
0245 %
0246 if fwd.source_ori == FIFF.FIFFV_MNE_FIXED_ORI || ...
0247         force_fixed == true
0248     nuse = 0;
0249     fwd.source_rr = zeros(fwd.nsource,3);
0250     fwd.source_nn = zeros(fwd.nsource,3);
0251     for k = 1:length(src)
0252         fwd.source_rr(nuse+1:nuse+src(k).nuse,:) = src(k).rr(src(k).vertno,:);
0253         fwd.source_nn(nuse+1:nuse+src(k).nuse,:) = src(k).nn(src(k).vertno,:);
0254         nuse = nuse + src(k).nuse;
0255     end
0256     %
0257     %   Modify the forward solution for fixed source orientations
0258     %
0259     if fwd.source_ori ~= FIFF.FIFFV_MNE_FIXED_ORI
0260         fprintf(1,'\tChanging to fixed-orientation forward solution...');
0261         fix_rot        = mne_block_diag(fwd.source_nn',1);
0262         fwd.sol.data   = fwd.sol.data*fix_rot;
0263         fwd.sol.ncol   = fwd.nsource;
0264         fwd.source_ori = FIFF.FIFFV_MNE_FIXED_ORI;
0265 
0266         if  ~isempty(fwd.sol_grad)
0267             fwd.sol_grad.data   = fwd.sol_grad.data*kron(fix_rot,eye(3));
0268             fwd.sol_grad.ncol   = 3*fwd.nsource;
0269         end
0270         fprintf(1,'[done]\n');
0271     end
0272 else
0273     if surf_ori
0274         %
0275         %   Rotate the local source coordinate systems
0276         %
0277         fprintf(1,'\tConverting to surface-based source orientations...');
0278         nuse = 0;
0279         pp   = 1;
0280         fwd.source_rr = zeros(fwd.nsource,3);
0281         for k = 1:length(src)
0282             fwd.source_rr(nuse+1:nuse+src(k).nuse,:) = src(k).rr(src(k).vertno,:);
0283             for p = 1:src(k).nuse
0284                 %
0285                 %  Project out the surface normal and compute SVD
0286                 %
0287                 nn = src(k).nn(src(k).vertno(p),:)';
0288                 [ U, S, V ]  = svd(eye(3,3) - nn*nn');
0289                 %
0290                 %  Make sure that ez is in the direction of nn
0291                 %
0292                 if nn'*U(:,3) < 0
0293                     U = -U;
0294                 end
0295                 fwd.source_nn(pp:pp+2,:) = U';
0296                 pp = pp + 3;
0297             end
0298             nuse = nuse + src(k).nuse;
0299         end
0300         surf_rot = mne_block_diag(fwd.source_nn',3);
0301         fwd.sol.data = fwd.sol.data*surf_rot;
0302         if  ~isempty(fwd.sol_grad)
0303             fwd.sol_grad.data = fwd.sol_grad.data*kron(surf_rot,eye(3));
0304         end
0305         fprintf(1,'[done]\n');
0306     else
0307         fprintf(1,'\tCartesian source orientations...');
0308         nuse = 0;
0309         fwd.source_rr = zeros(fwd.nsource,3);
0310         for k = 1:length(src)
0311             fwd.source_rr(nuse+1:nuse+src(k).nuse,:) = src(k).rr(src(k).vertno,:);
0312             nuse = nuse + src(k).nuse;
0313         end
0314         fwd.source_nn = kron(ones(fwd.nsource,1),eye(3,3));
0315         fprintf(1,'[done]\n');
0316     end
0317 end
0318 %
0319 % Add channel information
0320 %
0321 fwd.chs = chs;
0322 %
0323 %   Do the channel selection
0324 %
0325 if ~isempty(include) || ~isempty(exclude)
0326     %
0327     %   First do the channels to be included
0328     %
0329     if isempty(include)
0330         pick = ones(1,fwd.nchan);
0331     else
0332         pick = zeros(1,fwd.nchan);
0333         for k = 1:size(include,2)
0334             c = strmatch(include{k},fwd.sol.row_names,'exact');
0335             for p = 1:length(c)
0336                 pick(c(p)) = 1;
0337             end
0338         end
0339     end
0340     %
0341     %   Then exclude what needs to be excluded
0342     %
0343     if ~isempty(exclude)
0344         for k = 1:length(exclude)
0345             c = strmatch(exclude{k},fwd.sol.row_names,'exact');
0346             for p = 1:length(c)
0347                 pick(c(p)) = 0;
0348             end
0349         end
0350     end
0351     %
0352     %   Do we have something?
0353     %
0354     nuse = sum(pick);
0355     if nuse == 0
0356         error(me,'Nothing remains after picking');
0357     end
0358     %
0359     %   Create the selector
0360     %
0361     sel = zeros(1,nuse);
0362     p = 0;
0363     for k = 1:fwd.nchan
0364         if pick(k) > 0
0365             p = p + 1;
0366             sel(p) = k;
0367         end
0368     end
0369     fprintf(1,'\t%d out of %d channels remain after picking\n',nuse,fwd.nchan);
0370     %
0371     %   Pick the correct rows of the forward operator
0372     %
0373     fwd.nchan         = nuse;
0374     fwd.sol.data      = fwd.sol.data(sel,:);
0375     fwd.sol.nrow      = nuse;
0376     fwd.sol.row_names = fwd.sol.row_names(sel);
0377 
0378     if ~isempty(fwd.sol_grad)
0379         fwd.sol_grad.data      = fwd.sol_grad.data(sel,:);
0380         fwd.sol_grad.nrow      = nuse;
0381         fwd.sol_grad.row_names = fwd.sol_grad.row_names(sel);
0382     end
0383 
0384     fwd.chs = fwd.chs(sel);
0385 end
0386 
0387 return;
0388 
0389     function [tag] = find_tag(nodes,findkind)
0390 
0391    for kk = 1:length(nodes)
0392       node = nodes(kk);
0393       for p = 1:node.nent
0394          if node.dir(p).kind == findkind
0395             tag = fiff_read_tag(fid,node.dir(p).pos);
0396             return;
0397          end
0398       end
0399       tag = [];
0400       return;
0401    end
0402    
0403    end
0404 
0405     function [one] = read_one(node)
0406         %
0407         %   Read all interesting stuff for one forward solution
0408         %
0409         if isempty(node)
0410             one = [];
0411             return;
0412         end
0413         tag = find_tag(node,FIFF.FIFF_MNE_SOURCE_ORIENTATION);
0414         if isempty(tag)
0415             fclose(fid);
0416             error(me,'Source orientation tag not found');
0417         end
0418         one.source_ori = tag.data;
0419         tag = find_tag(node,FIFF.FIFF_MNE_COORD_FRAME);
0420         if isempty(tag)
0421             fclose(fid);
0422             error(me,'Coordinate frame tag not found');
0423         end
0424         one.coord_frame = tag.data;
0425         tag = find_tag(node,FIFF.FIFF_MNE_SOURCE_SPACE_NPOINTS);
0426         if isempty(tag)
0427             fclose(fid);
0428             error(me,'Number of sources not found');
0429         end
0430         one.nsource = tag.data;
0431         tag = find_tag(node,FIFF.FIFF_NCHAN);
0432         if isempty(tag)
0433             fclose(fid);
0434             error(me,'Number of channels not found');
0435         end
0436         one.nchan = tag.data;
0437         try
0438             one.sol = mne_transpose_named_matrix(fiff_read_named_matrix(fid,node,FIFF.FIFF_MNE_FORWARD_SOLUTION));
0439         catch
0440             fclose(fid);
0441             error(me,'Forward solution data not found (%s)',mne_omit_first_line(lasterr));
0442         end
0443         try
0444             one.sol_grad = mne_transpose_named_matrix(fiff_read_named_matrix(fid,node,FIFF.FIFF_MNE_FORWARD_SOLUTION_GRAD));
0445         catch
0446             one.sol_grad = [];
0447         end
0448         if size(one.sol.data,1) ~= one.nchan || ...
0449                 (size(one.sol.data,2) ~= one.nsource && size(one.sol.data,2) ~= 3*one.nsource)
0450             fclose(fid);
0451             error(me,'Forward solution matrix has wrong dimensions');
0452         end
0453         if ~isempty(one.sol_grad)
0454             if size(one.sol_grad.data,1) ~= one.nchan || ...
0455                     (size(one.sol_grad.data,2) ~= 3*one.nsource && size(one.sol_grad.data,2) ~= 3*3*one.nsource)
0456                 fclose(fid);
0457                 error(me,'Forward solution gradient matrix has wrong dimensions');
0458             end
0459         end
0460         return;
0461     end
0462 
0463 end