brainstormresults2freesurfer

 brainstormresults2freesurfer - script to extract brainstorm inverse to freesurfer

0001 % brainstormresults2freesurfer - script to extract brainstorm inverse to freesurfer
0002 
0003 % Loop over all controls and patients
0004 
0005 % Search ImageGridTime for selected time point
0006 % Extract Cdata from ImageGridAmp at selected time (column)
0007 
0008 %----------------------------------------------
0009 % Define component timing
0010 
0011 sa=[100 % first positive peak P100 (OT)
0012     150 % positive & negative peak P150 (SF) / N150 (OT)
0013     180 % negative N180 @ PT,AT and SP,SC
0014     250 % positive P250 (SP)
0015     400 ]; % P400 (SC,SP)
0016 
0017 wm=[ 80
0018     90 
0019     150
0020     260
0021     300
0022     400
0023     550 ];
0024 
0025 ea=[ 80 % N80  @ SF        , P80  @ OT/PT
0026     145   % N150 @ OT/PT     , P150 @ SPF (SC)
0027     240   %                  , P240 @ SPF/IPF & OT/SP
0028     300   % N300 @ L_PT/OT/IP,
0029     350   % Series: P350 @ Frontal (IPF), P500 @ SP
0030     450   % Series: P350 @ Frontal (IPF), P500 @ SP
0031     550 ];% Series: P350 @ Frontal (IPF), P500 @ SP; Also N550 @ left frontal
0032 
0033 %----------------------------------------------
0034 % Define groups
0035 groups = {'c','p'};
0036 
0037 %----------------------------------------------
0038 % Define subjects
0039 subs = [1:10];
0040 
0041 cond = {'ea','sa','wm'};
0042 
0043 for g=1:length(groups),
0044     for s=1:length(subs),
0045         
0046         sub = sprintf('%s%02d',groups{g},subs(s));
0047         
0048         %----------------------------------------------------
0049         % Load cortical tesselation from subjecttess.mat
0050         
0051         subtessfile = ['E:\matlab\brainstorm_v1\subjects\',sub,filesep,sub,'_subjecttess.mat'];
0052         
0053         if exist(subtessfile) ~= 2,
0054             % skip this group/subject, as no data available
0055             fprintf('\nskipping subject: %s\n\n',sub);
0056             continue;
0057         else
0058             p.mesh.path = ['E:\matlab\brainstorm_v1\subjects\',sub,filesep];
0059             p.mesh.file = [sub,'_subjecttess.mat'];
0060             p.mesh.type = 'BrainStorm';
0061            [p] = mesh_open(p);
0062             
0063             p.mesh.current = mesh_check(p,'cortex');
0064             
0065             % effectively remove all extraneous surfaces for mesh_write_freesurfer below
0066             indices = strcmp(p.mesh.data.meshtype,'cortex');
0067             exclude = find(indices == 0);
0068             p.mesh.data.meshtype(exclude) = [];
0069             
0070         end
0071         
0072         
0073         %----------------------------------------------------
0074         % define cortical source activity matrix file names
0075         
0076         %c01_oac_volts_data_results_image_svd12.mat
0077         %c01_ouc_volts_data_results_image_svd12.mat
0078         oucfile = ['E:\matlab\brainstorm_v1\studies\',sub,filesep,sub,'_ouc_volts_data_results_image_svd12.mat'];
0079         oacfile = ['E:\matlab\brainstorm_v1\studies\',sub,filesep,sub,'_oac_volts_data_results_image_svd12.mat'];
0080         oatfile = ['E:\matlab\brainstorm_v1\studies\',sub,filesep,sub,'_oat_volts_data_results_image_svd12.mat'];
0081         tacfile = ['E:\matlab\brainstorm_v1\studies\',sub,filesep,sub,'_tac_volts_data_results_image_svd12.mat'];
0082         
0083         if exist(oucfile) ~= 2, error(['no ',oucfile]); end;
0084         if exist(oacfile) ~= 2, error(['no ',oacfile]); end;
0085         if exist(oatfile) ~= 2, error(['no ',oatfile]); end;
0086         if exist(tacfile) ~= 2, error(['no ',tacfile]); end;
0087         
0088         %----------------------------------------------------
0089         % Save Cdata to freesurfer curvature files
0090         
0091         % D:\freesurfer\subjects\ptsdpet-c01\surf
0092         p.mesh.path = sprintf('D:\\freesurfer\\subjects\\ptsdpet-%s\\surf\\',sub);
0093         
0094         p.mesh.type = 'FS_curv'; % output freesurfer binary curvature files
0095         
0096         % cortical values are scaled into nano A.m (dipole moment)
0097         scale = 10^9;
0098         
0099         
0100         for c = 1:length(cond),
0101             
0102             
0103             switch cond{c},
0104                 
0105             case 'sa',
0106                 
0107                 %---------------------------------------------------------------------
0108                 % SA condition
0109                 
0110                 vert = size(p.mesh.data.vertices{p.mesh.current},1);
0111                 times = length(sa);
0112                 
0113                 Dif.oac = zeros(vert,times);
0114                 Dif.ouc = Dif.oac;
0115                 
0116                 load(oacfile,'ImageGridAmp','ImageGridTime'); % Time is in sec here
0117                 for t=1:length(sa),
0118                     
0119                     % find column of ImageGridAmp corresponding to wm(t)
0120                     col = find(ImageGridTime == sa(t)/1000);
0121                     if isempty(col), error('Cannot find time column in ImageGridTime\n'); end
0122                     
0123                     p.mesh.data.Cdata{p.mesh.current} = ImageGridAmp(:,col) .* scale;
0124                     Dif.oac(:,t) = p.mesh.data.Cdata{p.mesh.current};
0125                     
0126                     p.mesh.file = sprintf('rh.sa_%04d_oac_float',sa(t));
0127                     mesh_write(p);
0128                     
0129                 end
0130                 
0131                 clear ImageGridAmp ImageGridTime; pack;
0132                 
0133                 load(oucfile,'ImageGridAmp','ImageGridTime'); % Time is in sec here
0134                 for t=1:length(sa),
0135                     
0136                     % find column of ImageGridAmp corresponding to sa(t)
0137                     col = find(ImageGridTime == sa(t)/1000);
0138                     if isempty(col), error('Cannot find time column in ImageGridTime\n'); end
0139                     
0140                     p.mesh.data.Cdata{p.mesh.current} = ImageGridAmp(:,col) .* scale;
0141                     Dif.ouc(:,t) = p.mesh.data.Cdata{p.mesh.current};
0142                     
0143                     p.mesh.file = sprintf('rh.sa_%04d_ouc_float',sa(t));
0144                     mesh_write(p);
0145                 end
0146                 
0147                 clear ImageGridAmp ImageGridTime; pack;
0148                 
0149                 for t=1:length(sa),
0150                     
0151                     p.mesh.data.Cdata{p.mesh.current} = Dif.oac(:,t) - Dif.ouc(:,t);
0152                     p.mesh.file = sprintf('rh.sa_%04d_dif_float',sa(t));
0153                     mesh_write(p);
0154                     
0155                 end
0156                 
0157                 clear Dif; pack;
0158                 
0159                 
0160                 
0161             case 'wm',
0162                 
0163                 %---------------------------------------------------------------------
0164                 % WM condition
0165                 
0166                 vert = size(p.mesh.data.vertices{p.mesh.current},1);
0167                 times = length(wm);
0168                 
0169                 Dif.tac = zeros(vert,times);
0170                 Dif.oac = Dif.tac;
0171                 
0172                 load(tacfile,'ImageGridAmp','ImageGridTime'); % Time is in sec here
0173                 for t=1:length(wm),
0174                     
0175                     % find column of ImageGridAmp corresponding to wm(t)
0176                     col = find(ImageGridTime == wm(t)/1000);
0177                     if isempty(col), error('Cannot find time column in ImageGridTime\n'); end
0178                     
0179                     p.mesh.data.Cdata{p.mesh.current} = ImageGridAmp(:,col) .* scale;
0180                     Dif.tac(:,t) = p.mesh.data.Cdata{p.mesh.current};
0181                     
0182                     p.mesh.file = sprintf('rh.wm_%04d_tac_float',wm(t));
0183                     mesh_write(p);
0184                     
0185                 end
0186                 
0187                 clear ImageGridAmp ImageGridTime; pack;
0188                 
0189                 load(oacfile,'ImageGridAmp','ImageGridTime'); % Time is in sec here
0190                 for t=1:length(wm),
0191                     
0192                     % find column of ImageGridAmp corresponding to sa(t)
0193                     col = find(ImageGridTime == wm(t)/1000);
0194                     if isempty(col), error('Cannot find time column in ImageGridTime\n'); end
0195                     
0196                     p.mesh.data.Cdata{p.mesh.current} = ImageGridAmp(:,col) .* scale;
0197                     Dif.oac(:,t) = p.mesh.data.Cdata{p.mesh.current};
0198                     
0199                     p.mesh.file = sprintf('rh.wm_%04d_oac_float',wm(t));
0200                     mesh_write(p);
0201                 end
0202                 
0203                 clear ImageGridAmp ImageGridTime; pack;
0204                 
0205                 for t=1:length(wm),
0206                     
0207                     p.mesh.data.Cdata{p.mesh.current} = Dif.tac(:,t) - Dif.oac(:,t);
0208                     p.mesh.file = sprintf('rh.wm_%04d_dif_float',wm(t));
0209                     mesh_write(p);
0210                     
0211                 end
0212                 
0213                 clear Dif; pack;
0214                 
0215                 
0216                 
0217             case 'ea',
0218                 
0219                 %---------------------------------------------------------------------
0220                 % EA condition
0221                 
0222                 vert = size(p.mesh.data.vertices{p.mesh.current},1);
0223                 times = length(ea);
0224                 
0225                 Dif.oat = zeros(vert,times);
0226                 Dif.oac = Dif.oat;
0227                 
0228                 load(oatfile,'ImageGridAmp','ImageGridTime'); % Time is in sec here
0229                 for t=1:length(ea),
0230                     
0231                     % find column of ImageGridAmp corresponding to sa(t)
0232                     col = find(ImageGridTime == ea(t)/1000);
0233                     if isempty(col), error('Cannot find time column in ImageGridTime\n'); end
0234                     
0235                     p.mesh.data.Cdata{p.mesh.current} = ImageGridAmp(:,col) .* scale;
0236                     Dif.oat(:,t) = p.mesh.data.Cdata{p.mesh.current};
0237                     
0238                     p.mesh.file = sprintf('rh.ea_%04d_oat_float',ea(t));
0239                     mesh_write(p);
0240                     
0241                 end
0242                 
0243                 clear ImageGridAmp ImageGridTime; pack;
0244                 
0245                 load(oacfile,'ImageGridAmp','ImageGridTime'); % Time is in sec here
0246                 for t=1:length(ea),
0247                     
0248                     % find column of ImageGridAmp corresponding to sa(t)
0249                     col = find(ImageGridTime == ea(t)/1000);
0250                     if isempty(col), error('Cannot find time column in ImageGridTime\n'); end
0251                     
0252                     p.mesh.data.Cdata{p.mesh.current} = ImageGridAmp(:,col) .* scale;
0253                     Dif.oac(:,t) = p.mesh.data.Cdata{p.mesh.current};
0254                     
0255                     p.mesh.file = sprintf('rh.ea_%04d_oac_float',ea(t));
0256                     mesh_write(p);
0257                 end
0258                 
0259                 clear ImageGridAmp ImageGridTime; pack;
0260                 
0261                 for t=1:length(ea),
0262                     
0263                     p.mesh.data.Cdata{p.mesh.current} = Dif.oat(:,t) - Dif.oac(:,t);
0264                     p.mesh.file = sprintf('rh.ea_%04d_dif_float',ea(t));
0265                     mesh_write(p);
0266                     
0267                 end
0268                 
0269                 clear Dif; pack;
0270             end
0271         end
0272     end
0273 end