mesh_fit_elec - find mesh vertices nearest to electrodes
[p] = mesh_fit_elec(p)
This function is in development. A transformation matrix,
with rotations and translations, is in development to
facilitate coregistration of elec vertices to the scalp
vertices. This needs to be refined by a minimisation
function for the difference between the electrode
locations and the nearest scalp vertices. The scalp
mesh should be sufficiently refined to provide unique
vertices for each electrode vertex.
At present, it works with the BrainStorm meshes and the
elec_124_cart.txt file in the eeg_example_data folder.
The BrainStorm meshes contain a 'scalp' mesh with vertices
in meters and the electrode file contains position
vectors in centimeters. The electrodes are converted
to meters by ELEC_LOAD. In the example data, they were
previously coregistered and fitted to the scalp mesh.
The return vertices of the scalp that are nearest to
the electrodes are returned in the following:
p.mesh.data.meshtype{p.mesh.current}
p.mesh.data.vertices{p.mesh.current}
p.mesh.data.faces{p.mesh.current}
where p.mesh.current is also returned.
See also: fiducial_coregister0001 function [p] = mesh_fit_elec(p) 0002 0003 % mesh_fit_elec - find mesh vertices nearest to electrodes 0004 % 0005 % [p] = mesh_fit_elec(p) 0006 % 0007 % This function is in development. A transformation matrix, 0008 % with rotations and translations, is in development to 0009 % facilitate coregistration of elec vertices to the scalp 0010 % vertices. This needs to be refined by a minimisation 0011 % function for the difference between the electrode 0012 % locations and the nearest scalp vertices. The scalp 0013 % mesh should be sufficiently refined to provide unique 0014 % vertices for each electrode vertex. 0015 % 0016 % At present, it works with the BrainStorm meshes and the 0017 % elec_124_cart.txt file in the eeg_example_data folder. 0018 % The BrainStorm meshes contain a 'scalp' mesh with vertices 0019 % in meters and the electrode file contains position 0020 % vectors in centimeters. The electrodes are converted 0021 % to meters by ELEC_LOAD. In the example data, they were 0022 % previously coregistered and fitted to the scalp mesh. 0023 % 0024 % The return vertices of the scalp that are nearest to 0025 % the electrodes are returned in the following: 0026 % 0027 % p.mesh.data.meshtype{p.mesh.current} 0028 % p.mesh.data.vertices{p.mesh.current} 0029 % p.mesh.data.faces{p.mesh.current} 0030 % 0031 % where p.mesh.current is also returned. 0032 % 0033 % See also: fiducial_coregister 0034 % 0035 0036 % $Revision: 1.1 $ $Date: 2004/11/12 01:32:35 $ 0037 0038 % Licence: GNU GPL, no express or implied warranties 0039 % History: 04/2002, Darren.Weber_at_radiology.ucsf.edu 0040 % 09/2002, Darren.Weber_at_radiology.ucsf.edu 0041 % - bug fixing location of scalp mesh vertices 0042 % that correspond to electrode vertices. 0043 % - Still working with example dataset, not 0044 % a genuine coregistration algorithm 0045 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0046 0047 % debugging on/off 0048 global DB; 0049 DB = 0; 0050 test = 0; % coregistration testing 0051 0052 0053 0054 fprintf('MESH_FIT_ELEC: In development, not yet a true coregistration\n'); 0055 0056 p.mesh.current = mesh_check(p,'scalp'); 0057 0058 if isempty(p.mesh.current), 0059 error('MESH_FIT_ELEC: No ''scalp'' mesh in p struct.\n'); 0060 end 0061 0062 % Extract scalp triangulation 0063 FV.vertices = p.mesh.data.vertices{p.mesh.current}; 0064 FV.faces = p.mesh.data.faces{p.mesh.current}; 0065 0066 if isempty(p.elec.data), 0067 error('MESH_FIT_ELEC: ''p.elec.data'' is empty.\n'); 0068 else 0069 x = p.elec.data.x; 0070 y = p.elec.data.y; 0071 z = p.elec.data.z; 0072 end 0073 0074 fprintf('MESH_FIT_ELEC: ''Coregistering'' Electrodes...'); 0075 tic 0076 0077 if test, 0078 % Electrode fiducial points 0079 En = p.elec.data.nasion; 0080 Er = p.elec.data.rpa; 0081 El = p.elec.data.lpa; 0082 reg.fiducials.head = [ En; Er; El ]; 0083 0084 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0085 % NEED AN INTERACTIVE PROCESS TO SELECT THESE 0086 % Mesh fiducial points 0087 if isempty(p.mri.fiducials), 0088 if ~isempty(p.mri.data), 0089 fprintf('\n\nMESH_FIT_ELEC: Select Fiducial Points in MRI...\n\n'); 0090 avw_view(p.mri.data); 0091 waitfor(gcf); 0092 p.mri.fiducials = evalin('base',p.mri.fiducials); 0093 else 0094 fprintf('\n\nMESH_FIT_ELEC: Select Fiducial Points in MRI...\n\n'); 0095 [p] = mri_open(p); 0096 avw_view(p.mri.data); 0097 waitfor(gcf); 0098 p.mri.fiducials = evalin('base',p.mri.fiducials); 0099 end 0100 end 0101 reg.fiducials.mri = p.mri.fiducials; 0102 0103 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0104 T = fiducial_coregister(reg.fiducials.head,reg.fiducials.mri); 0105 0106 E = [p.elec.data.x p.elec.data.y p.elec.data.z]; 0107 E2MRI = [E,ones(size(E,1),1)] * T; 0108 end 0109 0110 0111 % perform least squares fit of elec to scalp??? 0112 %options = optimset('fminsearch'); 0113 %[sumd,fval,exitflag,output] = fminsearch('mesh_fit_elec_optim',... 0114 % 0, options, vertices, [x,y,z]); 0115 %fprintf('\n%s%f\n', 'Iterations = ', output.iterations); 0116 %fprintf('\n%s%d\n', 'Exit = ', exitflag); 0117 0118 0119 0120 0121 % Try to adjust the height of the electrodes to 0122 % those of the scalp mesh - really need a full 0123 % translation/rotation/scaling least squares solution 0124 elecmax = max(z); 0125 vertmax = max(FV.vertices(:,3)); 0126 z = z + (vertmax - elecmax); 0127 0128 0129 if DB, 0130 figure; plot3(x,y,z,'ro'); hold on; 0131 patch('vertices',FV.vertices,'faces',FV.faces,... 0132 'FaceColor',[0.0 0.0 0.6],'Edgecolor',[.6 .6 .6]); 0133 set(gca,'Projection','perspective'); 0134 set(gca,'DataAspectRatio',[1 1 1]); 0135 axis off tight vis3d 0136 end 0137 0138 0139 % Find nearest scalp vertices to electrode positions. 0140 % In the process, reorder the scalp mesh vertices 0141 % and correct the scalp mesh faces so that all scalp 0142 % vertices that lie near the electrodes 0143 % are numbered 1-N, in the same order as the 0144 % electrodes. This is done here so that the 0145 % mesh_laplacian_interp function will work 0146 % correctly. 0147 0148 % This mess might be better handled by adding 0149 % electrode vertices into the scalp mesh, rather 0150 % than finding the scalp mesh vertices nearest 0151 % to those of the electrodes. However, this would 0152 % require retesselation of the scalp mesh, which 0153 % is a tricky process 0154 0155 [FV,k] = mesh_nearest_reorder(FV,[x,y,z]); 0156 0157 0158 % Assign the adjusted triangulation back 0159 % into the original structure 0160 p.mesh.data.vertices{p.mesh.current} = FV.vertices; 0161 p.mesh.data.faces{p.mesh.current} = FV.faces; 0162 0163 0164 p.mesh.data.elecindex{p.mesh.current} = k; 0165 0166 x = FV.vertices(k,1); 0167 y = FV.vertices(k,2); 0168 z = FV.vertices(k,3); 0169 0170 % A quick triangulation of the scalp vertices 0171 % that lie nearest to the electrode vertices 0172 Efaces = convhulln(FV.vertices(k,:)); 0173 0174 % Now store the vertices 0175 [p.mesh.current,meshExists] = mesh_check(p,'elec'); 0176 0177 p.mesh.data.meshtype{p.mesh.current} = 'elec'; 0178 p.mesh.data.vertices{p.mesh.current} = [x,y,z]; 0179 p.mesh.data.faces{p.mesh.current} = Efaces; 0180 0181 0182 t = toc; fprintf('done (%6.2f sec).\n',t); 0183 0184 return 0185 0186 0187 0188 0189 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0190 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0191 function [FV,k] = mesh_nearest_reorder(FV,elecverts) 0192 0193 % Find the indices of the nearest vertices in 0194 % FV to those in vertices 0195 k = dsearchn(FV.vertices,elecverts); 0196 0197 if size(k,1) ~= size(elecverts,1), 0198 msg = sprintf('Sorry, duplicate vertices found in electrode/mesh coregistration'); 0199 error(msg); 0200 end 0201 0202 0203 % Debug visualization 0204 global DB; 0205 if DB, 0206 figure; hold on; 0207 patch('vertices',FV.vertices,'faces',FV.faces,... 0208 'FaceColor',[0.8 0.8 0.8],'Edgecolor',[.6 .6 .6]); 0209 set(gca,'Projection','perspective'); 0210 set(gca,'DataAspectRatio',[1 1 1]); 0211 axis off tight vis3d 0212 end 0213 0214 0215 % OK, now reorder the triangulation in FV 0216 % so that all the mesh vertices are ordered 0217 % 1:N 0218 0219 Nelec = size(elecverts,1); 0220 0221 for Eindex = 1:Nelec, 0222 0223 % It's important in this process that the 0224 % mesh indices are found over again, in case 0225 % they are modified for electrode Eindex-X 0226 elecpos = elecverts(Eindex,:); 0227 Kindex = dsearchn(FV.vertices,elecpos); 0228 0229 if Eindex == Kindex, 0230 continue; % great, do nothing 0231 end 0232 0233 % Get vertex coordinates at Eindex 0234 Evertex = FV.vertices(Eindex,:); 0235 % Get vertex coordinates at Kindex 0236 Kvertex = FV.vertices(Kindex,:); 0237 0238 if DB, 0239 fprintf('Fitting Electrode: %3d\n',Eindex); 0240 H1 = plot3(Evertex(1),Evertex(2),Evertex(3),'ro'); 0241 H2 = plot3(elecpos(1),elecpos(2),elecpos(3),'rd'); 0242 H3 = plot3(Kvertex(1),Kvertex(2),Kvertex(3),'bo'); 0243 pause 0244 delete(H1) 0245 delete(H2) 0246 delete(H3) 0247 end 0248 0249 % Now swap them 0250 FV.vertices(Kindex,:) = Evertex; 0251 FV.vertices(Eindex,:) = Kvertex; 0252 0253 % Swapping the vertices requires 0254 % rearranging the faces 0255 0256 % find all faces that contain Vindex 0257 EfaceIndices = find(FV.faces == Eindex); 0258 % find all faces that contain Kindex 0259 KfaceIndices = find(FV.faces == Kindex); 0260 0261 % Where faces refer to Kvertex, make them 0262 % refer to Vvertex and vice versa 0263 FV.faces(KfaceIndices) = Eindex; 0264 FV.faces(EfaceIndices) = Kindex; 0265 0266 % Debug testing.... 0267 if DB, 0268 % Get vertex coordinates at Vindex 0269 Evertex = FV.vertices(Eindex,:); 0270 % Get vertex coordinates at Kindex 0271 Kvertex = FV.vertices(Kindex,:); 0272 H1 = plot3(Evertex(1),Evertex(2),Evertex(3),'ro'); 0273 H2 = plot3(elecpos(1),elecpos(2),elecpos(3),'rd'); 0274 H3 = plot3(Kvertex(1),Kvertex(2),Kvertex(3),'bo'); 0275 pause 0276 0277 delete(H1) 0278 delete(H2) 0279 delete(H3) 0280 end 0281 0282 end 0283 0284 % Now the nearest vertices in FV to vertices 0285 % should be numbered 1:Nelec 0286 k = dsearchn(FV.vertices,elecverts); 0287 0288 if ~isequal(k,[1:Nelec]'), 0289 msg = sprintf('MESH_FIT_ELEC: Mesh vertices not in order of electrodes!'); 0290 warning(msg); 0291 end 0292 0293 return