Description of mesh_bem_shells

0001 
0002 % mesh_bem_shells_script - create bem shells for N subjects
0003 %
0004 % Assumes preprocessing of MRI with
0005 % FreeSurfer                http://surfer.nmr.mgh.harvard.edu/
0006 % FSL tools (BET & FAST)    http://fmrib.ox.ac.uk/
0007 
0008 %
0009 % Licence:  GNU GPL, no implied or express warranties
0010 % History:  08/2002, Darren.Weber_at_radiology.ucsf.edu
0011 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0012 
0013 
0014 
0015 S = 'c01';
0016 
0017 % These options are control flow points below
0018 meshplot    = 0;
0019 elecplot    = 0;
0020 getfid      = 0;
0021 coregister  = 0;
0022 
0023 % Create and initial correction of BEM
0024 % meshes from MRI and FreeSurfer data
0025 create      = 0;
0026 % Refinement of BEM correction process?
0027 correct     = 0;
0028 
0029 % default output is brainstorm format,
0030 % but these can be output also:
0031 % output BEM in EMSE format?
0032 emse        = 0;
0033 % output BEM in FreeSurfer format?
0034 freesurfer  = 1;
0035 
0036 
0037 
0038 % Fiducial points in RAS volume, obtained using avw_view (in meters)
0039 %
0040 %                  Nasion                 Right                Left
0041 %
0042 mriFID.sub{ 1} = 'c01';
0043 mriFID.xyz{ 1} = [  0.005  0.097 -0.022;  0.071  0.023 -0.052; -0.072  0.016 -0.050 ];
0044 mriFID.sub{ 2} = 'c02';
0045 mriFID.xyz{ 2} = [  0.007  0.088 -0.007;  0.081  0.015 -0.046; -0.070  0.012 -0.051 ];
0046 mriFID.sub{ 3} = 'c03';
0047 mriFID.xyz{ 3} = [ -0.002  0.098 -0.008;  0.055  0.040 -0.038; -0.095  0.030 -0.046 ];
0048 mriFID.sub{ 4} = 'c04';
0049 mriFID.xyz{ 4} = [  0.000  0.092  0.000;  0.070  0.024 -0.036; -0.076  0.031 -0.036 ];
0050 mriFID.sub{ 5} = 'c05';
0051 mriFID.xyz{ 5} = [  0.010  0.094 -0.010;  0.075  0.024 -0.048; -0.072  0.024 -0.052 ];
0052 mriFID.sub{ 6} = 'c06';
0053 mriFID.xyz{ 6} = [ -0.002  0.088  0.000;  0.077  0.008 -0.048; -0.075 -0.008 -0.048 ];
0054 mriFID.sub{ 7} = 'c07';
0055 mriFID.xyz{ 7} = [ -0.001  0.097 -0.018;  0.078  0.008 -0.055; -0.076  0.008 -0.055 ];
0056 mriFID.sub{ 8} = 'c08';
0057 mriFID.xyz{ 8} = [  0.000  0.096 -0.019;  0.087  0.005 -0.039; -0.065  0.004 -0.038 ];
0058 mriFID.sub{ 9} = 'c09';
0059 mriFID.xyz{ 9} = [  0.010  0.086 -0.020;  0.100 -0.020 -0.080; -0.060  0.020 -0.058 ];
0060 mriFID.sub{10} = 'c10';
0061 mriFID.xyz{10} = [  0.000  0.091 -0.018;  0.076  0.011 -0.052; -0.068  0.006 -0.052 ];
0062 
0063 mriFID.sub{11} = 'p02';
0064 mriFID.xyz{11} = [  0.002  0.091  0.016;  0.080  0.025 -0.036; -0.080  0.015 -0.043 ];
0065 mriFID.sub{12} = 'p04';
0066 mriFID.xyz{12} = [ -0.001  0.086 -0.009;  0.086  0.002 -0.062; -0.070  0.003 -0.065 ];
0067 mriFID.sub{13} = 'p05';
0068 mriFID.xyz{13} = [ -0.002  0.094 -0.001;  0.068  0.025 -0.040; -0.071  0.005 -0.042 ];
0069 mriFID.sub{14} = 'p06';
0070 mriFID.xyz{14} = [ -0.002  0.084  0.000;  0.082  0.013 -0.050; -0.066  0.013 -0.052 ];
0071 mriFID.sub{15} = 'p07';
0072 mriFID.xyz{15} = [  0.001  0.092  0.015;  0.080  0.003 -0.033; -0.070  0.004 -0.032 ];
0073 mriFID.sub{16} = 'p08';
0074 mriFID.xyz{16} = [ -0.003  0.095 -0.002;  0.070  0.018 -0.040; -0.074  0.022 -0.035 ];
0075 mriFID.sub{17} = 'p09';
0076 mriFID.xyz{17} = [ -0.002  0.100  0.004;  0.100  0.002 -0.028; -0.050  0.000 -0.036 ];
0077 
0078 
0079 
0080 
0081 
0082 if meshplot & ~correct,
0083     
0084     patch('vertices',p.mesh.data.vertices{4},'faces',p.mesh.data.faces{4},...
0085           'FaceColor',[1 0 0],'Edgecolor','none','FaceAlpha',.3);
0086     
0087     daspect([1 1 1]); axis tight; hold on
0088     
0089     patch('vertices',p.mesh.data.vertices{3},'faces',p.mesh.data.faces{3},...
0090           'FaceColor',[0 1 0],'Edgecolor','none','FaceAlpha',.5);
0091     
0092     patch('vertices',p.mesh.data.vertices{2},'faces',p.mesh.data.faces{2},...
0093           'FaceColor',[0 0 1],'Edgecolor','none','FaceAlpha',1.0);
0094     
0095     material dull
0096     camlight headlight
0097     
0098     mouse_rotate
0099     
0100     return
0101     
0102     % Check that vertices for scalp, oskull, iskull are coincident
0103     v = 1;
0104     
0105     vert = p.mesh.data.vertices{4};
0106     x = vert(v,1);
0107     y = vert(v,2);
0108     z = vert(v,3);
0109     plot3(x,y,z,'ro')
0110     
0111     vert = p.mesh.data.vertices{3};
0112     x = vert(v,1);
0113     y = vert(v,2);
0114     z = vert(v,3);
0115     plot3(x,y,z,'go')
0116     
0117     vert = p.mesh.data.vertices{2};
0118     x = vert(v,1);
0119     y = vert(v,2);
0120     z = vert(v,3);
0121     plot3(x,y,z,'bo')
0122     
0123     clear vert v x y z
0124     
0125     return
0126 end
0127 
0128 
0129 if elecplot,
0130     
0131     
0132     patch('vertices',p.mesh.data.vertices{4},'faces',p.mesh.data.faces{4},...
0133           'FaceColor',[1 0 0],'Edgecolor','none','FaceAlpha',.6);
0134     
0135     lighting phong
0136     material dull
0137     camlight headlight
0138     
0139     hold on
0140     
0141     x = p.elec.data.x;
0142     y = p.elec.data.y;
0143     z = p.elec.data.z;
0144     plot3(x,y,z,'bo')
0145     
0146 %     x = E2MRI(:,1);
0147 %     y = E2MRI(:,2);
0148 %     z = E2MRI(:,3);
0149 %     plot3(x,y,z,'bo')
0150     
0151     daspect([1 1 1]); axis tight; hold on
0152     
0153     mouse_rotate
0154     
0155     return
0156     
0157 end
0158 
0159 
0160 % 'c01', 'c02', 'c03', 'c04', 'c05', 'c06', 'c07', 'c08', 'c09', 'c10', 'p02', 'p04', 'p05', 'p06', 'p07', 'p08', 'p09'
0161 
0162 for sub = { S },
0163     
0164     % Image files to process
0165     
0166     % These images were created with the freesurfer mri_convert command, eg:
0167     % mri_convert -oid 1 0 0 -ojd 0 1 0 -okd 0 0 1 orig mri\analyze\c01_orig_axial_ras.img
0168     % These Analyze files were then processed with FSL tools to find the skull.
0169     IMG.path   = sprintf('\\\\POTZII\\data\\freesurfer\\subjects\\ptsdpet-%s\\mri\\analyze\\',char(sub));
0170     IMG.scalp  = sprintf('%s_orig_axial_ras',char(sub));
0171     IMG.oskull = sprintf('%s_orig_axial_ras_skull',char(sub));
0172     IMG.iskull = sprintf('%s_orig_axial_ras_bet',char(sub));
0173     IMG.intensity.scalp  = 100;
0174     IMG.intensity.oskull = 0.05;
0175     IMG.intensity.iskull = 150;
0176     IMG.tolerance.scalp  = 20;
0177     IMG.tolerance.oskull = 0.01;
0178     IMG.tolerance.iskull = 80;
0179     
0180     % Freesurfer surface to process
0181     % This freesurfer surface contains the whole brain surface.
0182     FS.path   = sprintf('\\\\POTZII\\data\\freesurfer\\subjects\\ptsdpet-%s\\surf\\',char(sub));
0183     FS.file   = 'rh.pial.asc';
0184     
0185     % Create the meshes and correct them (use pial cortex for this)
0186     if create,
0187        [p] = mesh_bem_shells(IMG,FS);
0188         p.mesh.path = FS.path;
0189         p.mesh.file = 'BS_subjecttess.mat';
0190         mesh_write(p);
0191        [p] = mesh_bem_correct(p);
0192         p.mesh.file = 'BS_corrected_subjecttess.mat';
0193         mesh_write(p);
0194         
0195         clear p
0196     end
0197     
0198     
0199     
0200     % Open already processed surfaces
0201    [p] = mesh_open; % init p struct
0202 %     p.mesh.path = FS.path;
0203 %     p.mesh.file = 'BS_subjecttess.mat';
0204 %     p.mesh.file = 'BS_corrected_subjecttess.mat';
0205     p.mesh.path = sprintf('\\\\POTZII\\data\\data_source\\%s\\meshes\\',char(sub));
0206     p.mesh.file = sprintf('%s_subjecttess.mat',char(sub));
0207    [p] = mesh_open(p);
0208     
0209     
0210 %     % LOAD SMOOTHWM AS CORTEX
0211 %     FS.file   = 'rh.smoothwm.asc';
0212 %     cd(FS.path)
0213 %     FS.surf = mesh_freesurfer2matlab(FS.file);
0214 %     Nsurf = 1;
0215 %     if isfield(FS.surf,'vertices'),
0216 %         p.mesh.data.meshtype{Nsurf} = 'cortex';
0217 %         p.mesh.data.vertices{Nsurf} = FS.surf.vertices;
0218 %         p.mesh.data.faces{Nsurf}    = FS.surf.faces;
0219 %         %patch('vertices',FS.surf.vertices,'faces',FS.surf.faces,'FaceColor',[.8 .0 .0],'Edgecolor','none');
0220 %     end
0221 %
0222 %     % Convert from mm to meters
0223 %     for m = 1:size(p.mesh.data.meshtype,2),
0224 %         p.mesh.data.vertices{m} = p.mesh.data.vertices{m} ./ 1000;
0225 %     end
0226 %     p.mesh.path = sprintf('\\\\POTZII\\data\\data_source\\%s\\meshes\\',char(sub));
0227 %     p.mesh.file = sprintf('%s_subjecttess.mat',char(sub));
0228 %     mesh_write(p);
0229     
0230     
0231     
0232     % Use manual revision to refine the dist struct and recorrect BEM
0233     if correct,
0234         dist.four2three = 5;    % scalp to outer skull
0235         dist.four2two   = 4;    % scalp to inner skull
0236         dist.three2two  = 3;    % outer skull to inner skull
0237         dist.two2one    = 1;    % inner skull to cortex
0238        [p] = mesh_bem_correct(p,[],dist,[],0);
0239     end
0240     
0241     
0242     
0243     
0244     % COREGISTER ELECTRODES TO BEM
0245     if coregister,
0246         
0247         % Use avw_view to update p.mri.fiducials in the base workspace
0248         % Then copy those values to the mriFID struct above and
0249         % adjust as necessary
0250         if getfid,
0251             cd(IMG.path)
0252             avw = avw_img_read(IMG.scalp);
0253             avw_view(avw);
0254             return
0255         end
0256         
0257         % Load the electrode data
0258         p.elec.path = sprintf('\\\\POTZII\\data\\data_source\\%s\\meshes\\',char(sub));
0259         p.elec.file = sprintf('%s_124fit.txt',char(sub));
0260         p.elec.plot = 0;
0261        [p] = elec_open(p);
0262         
0263         % Get surface fiducials from mriFID struct above
0264         Nfid = strmatch(sub,mriFID.sub);
0265         
0266         p.mri.fiducials = mriFID.xyz{Nfid}; % 3x3, nasion, right, left fiducial points in rows
0267         % Get electrode fiducials
0268         Efid = [p.elec.data.nasion; p.elec.data.rpa; p.elec.data.lpa];
0269         % Calculate coregistration transform
0270         T = fiducial_coregister(Efid,p.mri.fiducials);
0271         
0272         [path,file,ext] = fileparts([p.elec.path p.elec.file]);
0273         file = fullfile(path,[file,'_reg.mat']);
0274         save(file,'T')
0275         
0276         % Transform all the electrode coordinates to MRI space
0277         E = [p.elec.data.x p.elec.data.y p.elec.data.z];
0278         E2MRI = [E,ones(size(E,1),1)] * T;
0279         
0280         p.elec.data.x = E2MRI(:,1);
0281         p.elec.data.y = E2MRI(:,2);
0282         p.elec.data.z = E2MRI(:,3);
0283         
0284         REF2MRI = [p.elec.data.ref 1] * T;
0285         p.elec.data.ref    = REF2MRI(:,1:3);
0286         
0287         ORIGIN2MRI = [p.elec.data.centroid 1] * T;
0288         p.elec.data.centroid = ORIGIN2MRI(:,1:3);
0289         
0290         Efid2MRI = [ Efid [ 1 1 1]' ] * T;
0291         p.elec.data.nasion = Efid2MRI(1,1:3);
0292         p.elec.data.rpa    = Efid2MRI(2,1:3);
0293         p.elec.data.lpa    = Efid2MRI(3,1:3);
0294         
0295         
0296         clear E REF2MRI ORIGIN2MRI Efid2MRI
0297         
0298         %elec_write_3dspace(p);
0299         %elec_write_brainstorm(p);
0300         
0301         p.elec.file = sprintf('%s_124fit_new.elp',char(sub));
0302         elec_write_emse(p);
0303         
0304     end
0305     
0306     % Output BEM surfaces as EMSE format
0307     if emse,
0308         p.mesh.file = sprintf('%s.wfr',char(sub));
0309         p.mesh.type = 'emse';
0310         mesh_write(p);
0311     end
0312     
0313     % Output BEM surfaces as freesurfer format
0314     if freesurfer,
0315         p.mesh.path = FS.path;
0316         p.mesh.file = 'rh.asc';
0317         p.mesh.type = 'freesurfer';
0318         p.mesh.data.meshtype{1} = ''; % do not write cortex
0319         mesh_write(p);
0320     end
0321     
0322 end
0323 
0324 return
mesh_bem_shells_script

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
