eeg_interp_scalp_mesh

 eeg_interp_scalp_mesh - Script to interpolate scalp potentials over mesh
 
 This is part of the development and testing of
 functions to implement:

 Oostendorp T, Oosterom A, & Huiskamp G (1989),
 Interpolation on a triangulated 3D surface.  
 Journal of Computational Physics, 80: 331-343.

0001 % eeg_interp_scalp_mesh - Script to interpolate scalp potentials over mesh
0002 %
0003 % This is part of the development and testing of
0004 % functions to implement:
0005 %
0006 % Oostendorp T, Oosterom A, & Huiskamp G (1989),
0007 % Interpolation on a triangulated 3D surface.
0008 % Journal of Computational Physics, 80: 331-343.
0009 %
0010 
0011 % $Revision: 1.1 $ $Date: 2004/11/12 01:32:33 $
0012 
0013 % Licence:  GNU GPL, no implied or express warranties
0014 % History:  04/2002, Darren.Weber_at_radiology.ucsf.edu
0015 %           - needs verification and conversion to function
0016 %
0017 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0018 
0019 
0020 
0021 clear all
0022 
0023 % run R. Oostendorp's script to load simulated data
0024 cd d:\matlab\cvs\eeg_toolbox\lapint
0025 %sphere_load;
0026 %cd ..
0027 
0028 % If this script has already been run and the data
0029 % saved, then just load the saved data and return
0030 if isequal(exist('eeg_interp_scalp_mesh_data.mat'),2),
0031     fprintf('\nLoading saved data.\n');
0032     load 'eeg_interp_scalp_mesh_data';
0033     return
0034 end
0035 
0036 
0037 p = eeg_toolbox_defaults('create');
0038 % Load electrode coordinates 'Cartesian'
0039 p = elec_open(p);
0040 % Load potential values at each electrode
0041 p = eeg_open(p);
0042 % Load a tesselation dataset (scalp, 'skull', cortex)
0043 p = mesh_open(p);
0044 
0045 % Find the nearest scalp vertices to each
0046 % electrode vertex
0047 p = mesh_plot(p);
0048 close all
0049 
0050 % Extract tesselations from p struct
0051 scalpvert = p.meshData.vertices{3};
0052 scalpface = p.meshData.faces{3};
0053 elecvert  = p.meshData.vertices{4};
0054 elecface  = p.meshData.faces{4};
0055 
0056 % Extract electric potential at time t for electrodes
0057 Velec = p.voltageData;
0058 
0059 clear p;
0060 
0061 % This command retains the shape of scalp, while reducing
0062 % the number of faces from ~4000 to 1000.  In doing so, it
0063 % reduces the resolution of the scalp mesh and we need to
0064 % recompute the nearest scalp vertices for each electrode
0065 % vertex below
0066 %[scalpface, scalpvert] = reducepatch(scalpface,scalpvert,1000);
0067 
0068 % find the scalp vertex indices for each electrode vertex.
0069 % In this example, this generates replicate indices when
0070 % the reducepatch command above is executed
0071 Scalpindex = dsearchn(scalpvert,elecvert)';
0072 
0073 % Calculate the Laplacian matrix
0074 lap = mesh_laplacian(scalpvert,scalpface);
0075 
0076 % Calculate interpolation matrix, based on Laplacian
0077 [Lint, keepindex, repindex] = mesh_laplacian_interp(lap,Scalpindex);
0078 
0079 
0080 % Interpolate potential to scalp
0081 time = 100;
0082 if isempty(repindex),
0083     Vscalp = Lint * Velec(time,:)';
0084 else
0085     Vscalp = Lint * Velec(time,keepindex)';
0086 end
0087 
0088 %save eeg_interp_scalp_mesh_data
0089 
0090 
0091 return
0092 
0093 % Plot scalp with voltages at all interpolated vertices
0094 
0095 fig = figure;
0096 Hp = patch('vertices',scalpvert,'faces',scalpface,'FaceVertexCdata',Vscalp,'facecolor','interp','edgecolor','none');
0097 axis off tight vis3d
0098 lighting gouraud
0099 
0100 map = eeg_colormap('Red/Blue/White');
0101 colormap(map)
0102 
0103 [az,el] = view;
0104 lit = lightangle(az,el);
0105 
0106 set(Hp,'AmbientStrength',.7);
0107 set(Hp,'DiffuseStrength',.3);
0108 set(Hp,'SpecularStrength',0);
0109 set(Hp,'SpecularColorReflectance',0);
0110 
0111 %colorbar
0112 
0113 % Animation of the timeseries & save movie
0114 set(gcf,'BackingStore','off');
0115 set(Hp,'EraseMode','normal');
0116 figure(fig);
0117 % Define movie region as figure size
0118 rect = get(fig,'Position');
0119 rect(1:2) = [0 0];
0120 
0121 clear M
0122 
0123 start = 120;
0124 finish = 160;
0125 for t=start:finish,
0126     Vscalp = Lint * Velec(t,:)';
0127     set(Hp,'FaceVertexCdata',Vscalp);
0128     drawnow;
0129     M(t-(start-1)) = getframe(gcf,rect);
0130 end
0131 % Play the movie in another figure
0132 playtimes = 1;
0133 fps = 15;
0134 figure; axis off; movie(gcf,M,playtimes,fps,rect);
0135 movie2avi(M,'scalp_interpolation.avi','quality',100);
0136 
0137 
0138 
0139 return