extractresults2freesurfer

 Loop over all controls and patients

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