mesh_bem_correct

 mesh_bem_correct - Correct intersections of surfaces
 
[p] = mesh_bem_correct(p, origin, dist, coin, cortex)
 
 p - the eeg_toolbox_defaults struct; in this case it
 is important that it contains p.mesh.data with four cells, 
 corresponding, in order, of cortex, inner skull,
 outer skull and scalp meshes.  If any of these surfaces
 are not available, the corresponding cell should be
 left empty, but allocated.  See mesh_bem_shells
 to obtain these surfaces.  The vertex coordinates 
 should be given in meters.
 
 origin - 1x3, assumed to be [0 0 0] in Cartesian coordinates.
 
 dist - the minimum distance, in mm, between each layer of
 the boundary element model (BEM).  The function converts
 these values from mm to meters.  At present, it assumes
 four layers in the BEM (modify the code otherwise).  For 
 example:
 
     dist.four2three = 6;    % scalp to outer skull
     dist.four2two   = 8;    % scalp to inner skull
     dist.three2two  = 5;    % outer skull to inner skull
     dist.two2one    = 2;    % inner skull to cortex
 
 coin  - this function assumes that each vertex in scalp and
 skull layers are coincident (which they are when created
 with mesh_bem_shells).  If these surface vertices are not
 coincident, set this input parameter to zero.  The inner 
 skull to cortex vertices are assumed not to be coincident.
 However, note, at present, this function will not work for
 surfaces that are not coincident (this option may be enabled
 in future development, hack the code otherwise).
 
 cortex - 1 to process cortex/inner skull, 0 to process other
 surfaces only.  As the cortex/inner skull is the most time
 consuming process, this option allows for much faster 
 processing of the other surfaces.  This is useful for
 refined correction with adjustments of the dist struct.
 
 There are no MRI intensity checks in this function, it
 simply adjusts the shells so there are no intersections.  The
 inner skull is first located outside the cortex (given cortex
 argument is 1). The inner skull is then static.  Next, the 
 scalp is located a given distance outside the inner skull
 (usually moved at the inferior parietal and occipito-temporal 
 regions, if at all).  The scalp is then static.  The outer
 skull is then located some distance inside the scalp and
 then located some distance outside the inner skull.  A final
 check and iterative refined movement ensures that the outer 
 skull is located between the scalp and the inner skull.  Note
 that the scalp to outer skull distance can be larger than
 the final placement of the outer skull, as the outer skull
 is finally relocated with respect to the inner skull.

0001 function [p] = mesh_bem_correct(p,origin,dist,coin,cortex)
0002 
0003 % mesh_bem_correct - Correct intersections of surfaces
0004 %
0005 %[p] = mesh_bem_correct(p, origin, dist, coin, cortex)
0006 %
0007 % p - the eeg_toolbox_defaults struct; in this case it
0008 % is important that it contains p.mesh.data with four cells,
0009 % corresponding, in order, of cortex, inner skull,
0010 % outer skull and scalp meshes.  If any of these surfaces
0011 % are not available, the corresponding cell should be
0012 % left empty, but allocated.  See mesh_bem_shells
0013 % to obtain these surfaces.  The vertex coordinates
0014 % should be given in meters.
0015 %
0016 % origin - 1x3, assumed to be [0 0 0] in Cartesian coordinates.
0017 %
0018 % dist - the minimum distance, in mm, between each layer of
0019 % the boundary element model (BEM).  The function converts
0020 % these values from mm to meters.  At present, it assumes
0021 % four layers in the BEM (modify the code otherwise).  For
0022 % example:
0023 %
0024 %     dist.four2three = 6;    % scalp to outer skull
0025 %     dist.four2two   = 8;    % scalp to inner skull
0026 %     dist.three2two  = 5;    % outer skull to inner skull
0027 %     dist.two2one    = 2;    % inner skull to cortex
0028 %
0029 % coin  - this function assumes that each vertex in scalp and
0030 % skull layers are coincident (which they are when created
0031 % with mesh_bem_shells).  If these surface vertices are not
0032 % coincident, set this input parameter to zero.  The inner
0033 % skull to cortex vertices are assumed not to be coincident.
0034 % However, note, at present, this function will not work for
0035 % surfaces that are not coincident (this option may be enabled
0036 % in future development, hack the code otherwise).
0037 %
0038 % cortex - 1 to process cortex/inner skull, 0 to process other
0039 % surfaces only.  As the cortex/inner skull is the most time
0040 % consuming process, this option allows for much faster
0041 % processing of the other surfaces.  This is useful for
0042 % refined correction with adjustments of the dist struct.
0043 %
0044 % There are no MRI intensity checks in this function, it
0045 % simply adjusts the shells so there are no intersections.  The
0046 % inner skull is first located outside the cortex (given cortex
0047 % argument is 1). The inner skull is then static.  Next, the
0048 % scalp is located a given distance outside the inner skull
0049 % (usually moved at the inferior parietal and occipito-temporal
0050 % regions, if at all).  The scalp is then static.  The outer
0051 % skull is then located some distance inside the scalp and
0052 % then located some distance outside the inner skull.  A final
0053 % check and iterative refined movement ensures that the outer
0054 % skull is located between the scalp and the inner skull.  Note
0055 % that the scalp to outer skull distance can be larger than
0056 % the final placement of the outer skull, as the outer skull
0057 % is finally relocated with respect to the inner skull.
0058 %
0059 
0060 % $Revision: 1.1 $ $Date: 2004/11/12 01:32:35 $
0061 
0062 % Licence:  GNU GPL, no implied or express warranties
0063 % History:  10/2002, Darren.Weber_at_radiology.ucsf.edu
0064 %
0065 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0066 
0067 fprintf('\nMESH_BEM_CORRECT: ...\n'); tic;
0068 
0069 
0070 %fprintf('...sorry, undergoing repairs!\n\n'); return
0071 
0072 
0073 
0074 if ~exist('origin','var'), origin = [0 0 0];
0075 elseif isempty(origin),    origin = [0 0 0];
0076 end
0077 xo = origin(1); yo = origin(2); zo = origin(3);
0078 
0079 if ~exist('dist','var'),
0080     dist.four2three = 6;    % scalp to outer skull
0081     dist.four2two   = 8;    % scalp to inner skull
0082     dist.three2two  = 5;    % outer skull to inner skull
0083     dist.two2one    = 2;    % inner skull to cortex
0084 elseif isempty(dist),
0085     dist.four2three = 6;    % scalp to outer skull
0086     dist.four2two   = 8;    % scalp to inner skull
0087     dist.three2two  = 5;    % outer skull to inner skull
0088     dist.two2one    = 2;    % inner skull to cortex
0089 end
0090 
0091 % convert dist from mm to meters
0092 dist.four2three = dist.four2three / 1000;
0093 dist.four2two   = dist.four2two   / 1000;
0094 dist.three2two  = dist.three2two  / 1000;
0095 dist.two2one    = dist.two2one    / 1000;
0096 
0097 
0098 if ~exist('coin','var'),
0099     coin = 1;
0100 elseif isempty(dist),
0101     coin = 1;
0102 end
0103 
0104 if ~exist('cortex','var'),
0105     % run correction for inner skull/cortex
0106     cortex = 1;
0107 elseif isempty(cortex),
0108     cortex = 1;
0109 end
0110 
0111 % Identify each layer of the BEM
0112 [scalpN,exists]  = mesh_check(p,'scalp');
0113 if ~exists, error('...No scalp in p struct');
0114 else,       fprintf('...found cell %d of p.mesh.data contains scalp\n',scalpN);
0115 end
0116 [oskullN,exists] = mesh_check(p,'outer_skull');
0117 if ~exists, error('...No outer_skull in p struct');
0118 else,       fprintf('...found cell %d of p.mesh.data contains outer_skull\n',oskullN);
0119 end
0120 [iskullN,exists] = mesh_check(p,'inner_skull');
0121 if ~exists, error('...No inner_skull in p struct');
0122 else,       fprintf('...found cell %d of p.mesh.data contains inner_skull\n',iskullN);
0123 end
0124 [cortexN,exists] = mesh_check(p,'cortex');
0125 if ~exists,
0126     [cortexN,exists] = mesh_check(p,'pial');
0127     if ~exists, error('...No cortex in p struct');
0128     else fprintf('...found cell %d of p.mesh.data contains cortex\n',cortexN);
0129     end
0130 else,    fprintf('...found cell %d of p.mesh.data contains cortex\n',cortexN);
0131 end
0132 
0133 
0134 % Find distances and direction cosines of all surface vertices
0135 fprintf('...calculating vertex radius and direction cosines: scalp\n');
0136 Xscalp = p.mesh.data.vertices{scalpN}(:,1) - xo;
0137 Yscalp = p.mesh.data.vertices{scalpN}(:,2) - yo;
0138 Zscalp = p.mesh.data.vertices{scalpN}(:,3) - zo;
0139 Dscalp = sqrt( (Xscalp).^2 + (Yscalp).^2 + (Zscalp).^2 );
0140 Lscalp = (Xscalp)./Dscalp; % cos alpha
0141 Mscalp = (Yscalp)./Dscalp; % cos beta
0142 Nscalp = (Zscalp)./Dscalp; % cos gamma
0143 
0144 fprintf('...calculating vertex radius and direction cosines: outer skull\n');
0145 Xoskull = p.mesh.data.vertices{oskullN}(:,1) - xo;
0146 Yoskull = p.mesh.data.vertices{oskullN}(:,2) - yo;
0147 Zoskull = p.mesh.data.vertices{oskullN}(:,3) - zo;
0148 Doskull = sqrt( (Xoskull).^2 + (Yoskull).^2 + (Zoskull).^2 );
0149 Loskull = (Xoskull)./Doskull; % cos alpha
0150 Moskull = (Yoskull)./Doskull; % cos beta
0151 Noskull = (Zoskull)./Doskull; % cos gamma
0152 
0153 fprintf('...calculating vertex radius and direction cosines: inner skull\n');
0154 Xiskull = p.mesh.data.vertices{iskullN}(:,1) - xo;
0155 Yiskull = p.mesh.data.vertices{iskullN}(:,2) - yo;
0156 Ziskull = p.mesh.data.vertices{iskullN}(:,3) - zo;
0157 Diskull = sqrt( (Xiskull).^2 + (Yiskull).^2 + (Ziskull).^2 );
0158 Liskull = (Xiskull)./Diskull; % cos alpha
0159 Miskull = (Yiskull)./Diskull; % cos beta
0160 Niskull = (Ziskull)./Diskull; % cos gamma
0161 
0162 % Note that the L,M,N direction cosines of the above
0163 % should be equal when these surface vertices are
0164 % coincident (given 1-2 decimal places).
0165 
0166 % check that max distances are less than 1 meter
0167 if max(Dscalp) > 1,  error('...maximum scalp radius > 1 meter.\n'); end
0168 if max(Doskull) > 1, error('...maximum outer skull radius > 1 meter.\n'); end
0169 if max(Diskull) > 1, error('...maximum inner skull radius > 1 meter.\n'); end
0170 
0171 
0172 CORRECTION_ADJUSTMENT = 0.0001;
0173 
0174 
0175 if cortex,
0176     
0177     iSkull.vertices = p.mesh.data.vertices{iskullN};
0178     iSkull.faces    = p.mesh.data.faces{iskullN};
0179     
0180     cortex.vertices = p.mesh.data.vertices{cortexN};
0181     cortex.faces    = p.mesh.data.faces{cortexN};
0182     
0183     fprintf('...processing cortex/inner skull (a lengthy process)...\n');
0184     fprintf('...tesselating cortex vertices for searching...'); tic;
0185     tri = delaunayn(cortex.vertices);
0186     t = toc; fprintf(' (%5.2f sec)\n',t);
0187     
0188     correct = 1;
0189     while correct,
0190         % Find cortex vertices closest to the inner skull vertex
0191         iSkullNvert = size(iSkull.vertices,1);
0192         fprintf('...checking for nearest cortex vertex to %d inner skull vertices...',iSkullNvert); tic;
0193         cortex_vertex_index = dsearchn(cortex.vertices,tri,iSkull.vertices);
0194         t = toc; fprintf(' (%5.2f sec)\n',t);
0195         
0196         Xcortex = cortex.vertices(cortex_vertex_index,1) - xo;
0197         Ycortex = cortex.vertices(cortex_vertex_index,2) - yo;
0198         Zcortex = cortex.vertices(cortex_vertex_index,3) - zo;
0199         Dcortex = sqrt( (Xcortex).^2 + (Ycortex).^2 + (Zcortex).^2 );
0200         if max(Dcortex) > 1, error('...maximum cortex radius > 1 meter.\n'); end
0201         
0202         % Don't calculate direction cosines for cortex, because these
0203         % vertices are not to be moved
0204         
0205         % ensure the inner skull vertices are outside the cortex
0206         dif = Diskull - Dcortex;
0207         correct = find(dif < dist.two2one);
0208         
0209         if correct,
0210             
0211             % plot the cortex and inner skull
0212             iSkullfig = figure;
0213             patch('vertices',cortex.vertices,'faces',cortex.faces,'facecolor',[.7 0  0],'edgecolor','none'); hold on
0214             patch('vertices',iSkull.vertices,'faces',iSkull.faces,'facecolor',[ 0 0 .7],'edgecolor','none','facealpha',.4);
0215             axis vis3d; view(3); axis tight; light
0216             drawnow
0217             
0218             vert = iSkull.vertices(correct,:);
0219             x = vert(:,1);
0220             y = vert(:,2);
0221             z = vert(:,3);
0222             plot3(x,y,z,'bo')
0223             
0224             pause(5)
0225             close(iSkullfig);
0226             
0227             fprintf('...correcting %4d inner skull vertices < cortex + %d mm\n',size(correct,1),dist.two2one * 1000);
0228             Diskull(correct) = Diskull(correct) + (dist.two2one - dif(correct)) + CORRECTION_ADJUSTMENT;
0229             Xiskull = (Liskull .* Diskull);
0230             Yiskull = (Miskull .* Diskull);
0231             Ziskull = (Niskull .* Diskull);
0232             iSkull.vertices = [ Xiskull Yiskull Ziskull ];
0233             iSkull = mesh_smooth(iSkull,origin);
0234             
0235         else
0236             fprintf('...all inner skull vertices > cortex + %d mm\n',dist.two2one * 1000);
0237         end
0238         % iterate until corrected, finding cortex vertices
0239         % each time, as the nearest cortex vertex can change
0240         if ishandle(iSkullfig), close(iSkullfig); end
0241     end
0242     
0243     p.mesh.data.vertices{iskullN} = iSkull.vertices;
0244     p.mesh.data.faces{iskullN}    = iSkull.faces;
0245     clear iSkull cortex;
0246     
0247 end
0248 
0249 % ensure the scalp vertices are outside the inner skull
0250 FV.vertices = p.mesh.data.vertices{scalpN};
0251 FV.faces    = p.mesh.data.faces{scalpN};
0252 
0253 correct = 1;
0254 while correct,
0255     dif = Dscalp - Diskull;
0256     correct = find(dif < dist.four2two);
0257     if correct,
0258         fprintf('...correcting %4d scalp vertices < inner skull + %d mm\n',size(correct,1),dist.four2two * 1000);
0259         Dscalp(correct) = Dscalp(correct) + (dist.four2two - dif(correct)) + CORRECTION_ADJUSTMENT;
0260         Xscalp = (Lscalp .* Dscalp);
0261         Yscalp = (Mscalp .* Dscalp);
0262         Zscalp = (Nscalp .* Dscalp);
0263         FV.vertices = [ Xscalp Yscalp Zscalp ];
0264         FV = mesh_smooth(FV,origin);
0265     else
0266         fprintf('...all scalp vertices > inner skull + %d mm\n',dist.four2two * 1000);
0267     end
0268 end
0269 p.mesh.data.vertices{scalpN} = FV.vertices;
0270 p.mesh.data.faces{scalpN}    = FV.faces;
0271 clear FV;
0272 
0273 % ensure the outer skull vertices are inside the scalp
0274 FV.vertices = p.mesh.data.vertices{oskullN};
0275 FV.faces    = p.mesh.data.faces{oskullN};
0276 
0277 correct = 1;
0278 while correct,
0279     dif = Dscalp - Doskull;
0280     correct = find(dif < dist.four2three);
0281     if correct,
0282         fprintf('...correcting %4d outer skull vertices > scalp - %d mm\n',size(correct,1),dist.four2three * 1000);
0283         Doskull(correct) = Doskull(correct) - (dist.four2three - dif(correct)) - CORRECTION_ADJUSTMENT;
0284         Xoskull = (Loskull .* Doskull);
0285         Yoskull = (Moskull .* Doskull);
0286         Zoskull = (Noskull .* Doskull);
0287         FV.vertices = [ Xoskull Yoskull Zoskull ];
0288         FV = mesh_smooth(FV,origin);
0289     else
0290         fprintf('...all outer skull vertices < scalp - %d mm\n',dist.four2three);
0291     end
0292 end
0293 
0294 % ensure the outer skull vertices are outside the inner skull
0295 correct = 1;
0296 while correct,
0297     dif = Doskull - Diskull;
0298     correct = find(dif < dist.three2two);
0299     if correct,
0300         fprintf('...correcting %4d outer skull vertices < inner skull + %d mm\n',size(correct,1),dist.three2two * 1000);
0301         Doskull(correct) = Doskull(correct) + (dist.three2two - dif(correct)) + CORRECTION_ADJUSTMENT;
0302         Xoskull = (Loskull .* Doskull);
0303         Yoskull = (Moskull .* Doskull);
0304         Zoskull = (Noskull .* Doskull);
0305         FV.vertices = [ Xoskull Yoskull Zoskull ];
0306         FV = mesh_smooth(FV,origin);
0307     else
0308         fprintf('...all outer skull vertices > inner skull + %d mm\n',dist.three2two * 1000);
0309     end
0310 end
0311 
0312 
0313 % finally, ensure the outer skull vertices are between the scalp and the inner skull
0314 measure = .002; %start with 2 mm
0315 pass = 0;
0316 move = [1 1];
0317 while find(move),
0318     
0319     pass = pass + 1;
0320     if pass > 1,
0321         % decrease the testing measure, to allow progressive
0322         % determination of outer skull within scalp/inner skull
0323         measure = measure - 0.00025;
0324     end
0325     
0326     % Confirm that oskull is at least 2mm inside scalp
0327     scalp2oskull = Dscalp - Doskull;
0328     correct = find(scalp2oskull < measure);
0329     if correct,
0330         fprintf('...correcting %4d outer skull vertices > scalp - %d mm\n',size(correct,1),measure * 1000);
0331         Doskull(correct) = Doskull(correct) - (measure - scalp2oskull(correct)) - CORRECTION_ADJUSTMENT;
0332         Xoskull = (Loskull .* Doskull);
0333         Yoskull = (Moskull .* Doskull);
0334         Zoskull = (Noskull .* Doskull);
0335         FV.vertices = [ Xoskull Yoskull Zoskull ];
0336         move(1) = 1;
0337     else,
0338         fprintf('...confirmed outer skull vertices < scalp - %d mm\n',measure * 1000);
0339         move(1) = 0;
0340     end
0341     
0342     oskull2iskull = Doskull - Diskull;
0343     correct = find(oskull2iskull < measure);
0344     if correct,
0345         fprintf('...correcting %4d outer skull vertices < inner skull + %d mm\n',size(correct,1),measure * 1000);
0346         Doskull(correct) = Doskull(correct) + (measure - oskull2iskull(correct)) + CORRECTION_ADJUSTMENT;
0347         Xoskull = (Loskull .* Doskull);
0348         Yoskull = (Moskull .* Doskull);
0349         Zoskull = (Noskull .* Doskull);
0350         FV.vertices = [ Xoskull Yoskull Zoskull ];
0351         move(2) = 1;
0352     else,
0353         fprintf('...confirmed outer skull vertices > inner skull + %d mm\n',measure * 1000);
0354         move(2) = 0;
0355     end
0356 end
0357 p.mesh.data.vertices{oskullN} = FV.vertices;
0358 p.mesh.data.faces{oskullN}    = FV.faces;
0359 clear FV;
0360 
0361 
0362 t=toc; fprintf('...done (%5.2f sec).\n\n',t);
0363 
0364 return
0365 
0366 
0367 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0368 % This code can be used to plot the vertices that are to be
0369 % corrected
0370 %
0371 %     v = correct;
0372 %     index = [scalpN oskullN];
0373 %
0374 %     patch('vertices',p.mesh.data.vertices{index(1)},'faces',p.mesh.data.faces{index(1)},...
0375 %           'FaceColor',[.9 .9 .9],'Edgecolor','none','FaceAlpha',.6);
0376 %     daspect([1 1 1]); axis tight; hold on
0377 %
0378 %     vert = p.mesh.data.vertices{index(1)};
0379 %     x = vert(v,1);
0380 %     y = vert(v,2);
0381 %     z = vert(v,3);
0382 %     plot3(x,y,z,'ro')
0383 %
0384 %     patch('vertices',p.mesh.data.vertices{index(2)},'faces',p.mesh.data.faces{index(2)},...
0385 %           'FaceColor',[.0 .6 .0],'Edgecolor',[.2 .2 .2],'FaceAlpha',.8);
0386 %
0387 %     vert = p.mesh.data.vertices{index(2)};
0388 %     x = vert(v,1);
0389 %     y = vert(v,2);
0390 %     z = vert(v,3);
0391 %     plot3(x,y,z,'bo')