mesh_scalp_interp

 mesh_scalp_interp - Interpolates scalar over scalp mesh from electrodes
 
 Usage: [p] = mesh_scalp_interp(p)
 
 p is the eeg_toolbox struct (see eeg_toolbox_defaults)
 
 This function returns the laplacian matrix of the scalp
 mesh, the interpolation matrix based on the min norm of
 the laplacian, and the voltage interpolated on the scalp.
 It also calls an electrode "coregistration" function and
 the results of that operation.
 
 The return matrices are all in p.mesh.data, including the 
 interpolated timeseries values in p.mesh.data.timeseries{N},
 where N is the index of the scalp mesh in this case (see 
 mesh_check to get N).
 
 This function is called by eeg_contours_engine.
 
 See also, mesh_fit_elec, mesh_laplacian, mesh_laplacian_interp
 for more information.

0001 function [p] = mesh_scalp_interp(p)
0002 
0003 % mesh_scalp_interp - Interpolates scalar over scalp mesh from electrodes
0004 %
0005 % Usage: [p] = mesh_scalp_interp(p)
0006 %
0007 % p is the eeg_toolbox struct (see eeg_toolbox_defaults)
0008 %
0009 % This function returns the laplacian matrix of the scalp
0010 % mesh, the interpolation matrix based on the min norm of
0011 % the laplacian, and the voltage interpolated on the scalp.
0012 % It also calls an electrode "coregistration" function and
0013 % the results of that operation.
0014 %
0015 % The return matrices are all in p.mesh.data, including the
0016 % interpolated timeseries values in p.mesh.data.timeseries{N},
0017 % where N is the index of the scalp mesh in this case (see
0018 % mesh_check to get N).
0019 %
0020 % This function is called by eeg_contours_engine.
0021 %
0022 % See also, mesh_fit_elec, mesh_laplacian, mesh_laplacian_interp
0023 % for more information.
0024 %
0025 
0026 % $Revision: 1.1 $ $Date: 2004/11/12 01:32:35 $
0027 
0028 % Licence:  Gnu GPL, no express or implied warranties
0029 % History:  09/02, Darren.Weber_at_radiology.ucsf.edu
0030 %                  extracted out of MESH_SCALP_INTERP, called by that
0031 %                  function.
0032 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0033 
0034 if ~exist('p','var'),
0035     fprintf('Setting default parameters.\n');
0036    [p] = eeg_toolbox_defaults('create');
0037 end
0038 
0039 % Calculate voltage at scalp mesh vertices
0040 if isempty(p.mesh.data),
0041     if ~isempty(p.mesh.file),
0042        [p] = mesh_open(p);
0043     else
0044         msg = sprintf('MESH_SCALP_INTERP: No meshes, load meshes first.');
0045         warning(msg);
0046     end
0047 end
0048 
0049 % get the scalp mesh
0050 [p.mesh.current,meshExists] = mesh_check(p,'scalp');
0051 
0052 if isempty(p.mesh.current),
0053     msg = sprintf('MESH_SCALP_INTERP: No ''scalp'' mesh, check p.mesh.data.meshtype.');
0054     error(msg);
0055 end
0056 
0057 % Find nearest vertices of scalp mesh to electrodes
0058 % This will also reorder the scalp mesh vertices/faces
0059 % so that the first 1:Nelec vertices are the vertices
0060 % in the scalp that lie closest to the electrodes.
0061 % Serious problems arise if there are any duplicate
0062 % scalp vertices for any electrode, but this routine
0063 % will continue regardless, at present.
0064 p = mesh_fit_elec(p);
0065 
0066 % Ensure scalp is the currently selected mesh
0067 elecN = mesh_check(p,'elec');
0068 scalpN = mesh_check(p,'scalp');
0069 [p.mesh.current,meshExists] = mesh_check(p,'scalp');
0070 
0071 % Get vertices and faces of scalp mesh, after they
0072 % have been rearraned by mesh_fit_elec
0073 scalpvert = p.mesh.data.vertices{scalpN};
0074 scalpface = p.mesh.data.faces{scalpN};
0075 
0076 % Get the scalp vertices and faces
0077 % that correspond to the electrodes
0078 elecvert  = p.mesh.data.vertices{elecN};
0079 elecface  = p.mesh.data.faces{elecN};
0080 % Get the index of each electrode vertex
0081 % into the scalp mesh, should be first
0082 % 1:Nelec indices of scalp vertices
0083 elecindex = p.mesh.data.elecindex{scalpN};
0084 
0085 fprintf('MESH_SCALP_INTERP: Calculating Scalp Interpolation\n');
0086 
0087 
0088 % Could identify the elements of the scalp below
0089 % the electrodes and discard them before calculating
0090 % the interpolation matrices, or use them in the
0091 % interpolation, but set their potential to zero
0092 % later (see comment below).
0093 
0094 
0095 % Calculate the Laplacian matrix of the scalp mesh (once off calculation)
0096 if isfield(p.mesh.data,'lapmat'),
0097     if isempty(p.mesh.data.lapmat{scalpN}),
0098         p.mesh.data.lapmat{scalpN} = mesh_laplacian(scalpvert,scalpface);
0099     end
0100 else
0101     p.mesh.data.lapmat{scalpN} = mesh_laplacian(scalpvert,scalpface);
0102 end
0103 
0104 
0105 % Calculate the electrode to scalp mesh interpolation matrix
0106 % using a Laplacian minimum norm constraint (once off calculation)
0107 if isfield(p.mesh.data,'lapint'),
0108     if isempty(p.mesh.data.lapint{scalpN}),
0109         [p.mesh.data.lapint{scalpN}, p.mesh.data.keepindex{scalpN}, p.mesh.data.repindex{scalpN}] = mesh_laplacian_interp(p.mesh.data.lapmat{scalpN}, elecindex);
0110     end
0111 else
0112     [p.mesh.data.lapint{scalpN}, p.mesh.data.keepindex{scalpN}, p.mesh.data.repindex{scalpN}] = mesh_laplacian_interp(p.mesh.data.lapmat{scalpN}, elecindex);
0113 end
0114 
0115 
0116 fprintf('MESH_SCALP_INTERP: Multiplying Voltage by Interpolation matrix...\n');
0117 tic;
0118 % Interpolate voltage to all scalp mesh vertices (once off calc)
0119 if isempty(p.mesh.data.repindex{scalpN}),
0120     p.mesh.data.Cdata{scalpN} = p.mesh.data.lapint{scalpN} * p.volt.data';
0121     %p.mesh.data.Cdata{scalpN} = p.volt.data * p.mesh.data.lapint{scalpN}';
0122 else
0123     p.mesh.data.Cdata{scalpN} = p.mesh.data.lapint{scalpN} * p.volt.data(:,p.mesh.data.keepindex{scalpN})';
0124     %p.mesh.data.Cdata{scalpN} = p.volt.data(:,p.mesh.data.keepindex{scalpN}) * p.mesh.data.lapint{scalpN}';
0125 end
0126 t = toc;
0127 fprintf('done (%6.2f sec).\n',t);
0128 
0129 
0130 % Now a quick test of the interpolation: the values at the first
0131 % Nelec vertices of the scalp mesh must be equal to the values
0132 % of the original voltage data
0133 TMP = p.mesh.data.Cdata{scalpN};
0134 Nelec = size(p.volt.data,2);
0135 if ~isequal(p.volt.data',TMP(1:Nelec,:)),
0136     msg = sprintf('MESH_SCALP_INTERP: Fatal Error in Mesh Voltage Interpolation, Reload all Data?\n');
0137     error(msg);
0138 end
0139 
0140 
0141 % At present, all scalp vertices are interpolated.  It might
0142 % be useful at this point to determine all vertices that are
0143 % "well below" the electrode array, say > 5cm from a nearest
0144 % electrode and then set their voltage to zero.
0145 
0146 
0147 
0148 return