timecourse

 timecourse.m

 Extracts time-course from an ROI

0001 % timecourse.m
0002 %
0003 % Extracts time-course from an ROI
0004 %
0005 
0006 global BCH; %- used as a flag to know if we are in batch mode or not.
0007 
0008 USE_INTERFACE = 1;
0009 USE_YMAD = 0;
0010 
0011 rootpath = '/home03/handy/imagery/';
0012 spmpath = 'analysis';
0013 sub_ids = ...
0014     {'02jun00CW', ...
0015      '03aug00GG', ...
0016      '10aug00EH', ...
0017      '11aug00WS', ...
0018      '14jun00OMG', ...
0019      '15jun00JE', ...
0020      '15jun00SN', ...
0021      '17jul00RC', ...
0022      '18jul00CS'}; 
0023 
0024 series_map.episodic = [ ...
0025     3 6; ...
0026     2 5; ...
0027     2 6; ...
0028     3 6; ...
0029     3 6; ...
0030     2 5; ...
0031     3 6; ...
0032     2 5; ...
0033     3 6];
0034 
0035 series_map.semantic = [ ...
0036     1 4; ...
0037     3 6; ...
0038     3 6; ...
0039     1 4; ...
0040     1 4; ...
0041     3 6; ...
0042     1 4; ...
0043     3 6; ...
0044     1 4];
0045 
0046 task_vols = [16:34 50:68 84:102 118:136];
0047 rest_vols = [1:15 35:49 69:83 103:117];
0048 
0049 nser = 6;
0050 nsub = length(sub_ids);
0051 
0052 
0053 if USE_INTERFACE
0054 
0055   %
0056   %  Deal with getting the cluster of interest from the image of interest.  This
0057   %  is the cluster that will be used to pull data out of individual subject files.
0058   %
0059 
0060   % from spm_regions.m
0061 
0062   Finter = spm_figure('GetWin','Interactive');
0063   Fgraph = spm_figure('GetWin','Graphics');
0064   set(Finter,'Name','VOI time-series extraction')
0065 
0066   %-Find nearest voxel [Euclidean distance] in point list with data saved
0067   % in Y.mad, and update GUI location
0068   %-----------------------------------------------------------------------
0069   if ~length(SPM.XYZmm)
0070     spm('alert!','No suprathreshold voxels!',mfilename,0);
0071     Y = []; xY = [];
0072     return
0073   elseif exist(fullfile(SPM.swd,'Y.mad')) ~= 2
0074     spm('alert!','No raw data saved with this analysis!',mfilename,0);
0075     Y = []; xY = [];
0076     return
0077   elseif ~any(SPM.QQ)
0078     spm('alert!','No raw data saved for any suprathreshold location!',...
0079     mfilename,0);
0080     Y = []; xY = [];
0081     return
0082   end
0083 
0084   % from spm_VOI.m and spm_regions.m
0085   xyzmm   = spm_results_ui('GetCoords');
0086 
0087   [xyzmm,i] = spm_XYZreg('NearestXYZ',xyzmm,SPM.XYZmm);
0088   spm_results_ui('SetCoords',xyzmm);
0089   A     = spm_clusters(SPM.XYZ);
0090   Q     = find(A == A(i));
0091   str   = sprintf('%0.0f voxel cluster',length(Q))
0092   D     = SPM.XYZ(:,Q);
0093   S     = length(Q);
0094   
0095   q       = find(SPM.QQ(Q));
0096 
0097   if any(SPM.QQ(Q)==0)
0098     spm('alert"',{...
0099         sprintf('Don''t have raw data for all %d suprathreshold',length(Q)),...
0100     sprintf('voxels.'), ...
0101     sprintf('Proceeding using the %d voxels that are.',length(q)),...
0102          },mfilename,sqrt(-1));
0103   end
0104 else
0105   D = global_cluster;
0106 end % if USE_INTERFACE
0107 
0108 %
0109 % Read the data from one of two sources.
0110 %   One is the Y.mad file
0111 %   The other is from raw data files
0112 
0113 %-Get (approximate) raw data y from Y.mad file
0114 %-NB: Data in Y.mad file is compressed, and therefore not fully accurate
0115 %-----------------------------------------------------------------------
0116 
0117 clear y
0118 
0119 if USE_YMAD
0120   y       = spm_extract(fullfile(SPM.swd,'Y.mad'),SPM.QQ(Q(q)));
0121 else
0122 
0123   for isub = 1:nsub
0124     clear tmp
0125     for iser = 1:nser
0126       switch sub_ids{isub}
0127        case '14jun00OMG'
0128     tmp{iser} = spm_get('Files', ... 
0129                 sprintf('%s%s/spm/functional/series%02d/', ...
0130                     rootpath, sub_ids{isub}, iser), ...
0131                 'sns*14jun000MG*.img');
0132        otherwise
0133     tmp{iser} = spm_get('Files', ... 
0134                 sprintf('%s%s/spm/functional/series%02d/', rootpath, sub_ids{isub}, iser), ... 
0135                 sprintf('sns*%s*.img', sub_ids{isub}));
0136       end
0137     end
0138     filelist{isub} = char(tmp);
0139 
0140 %    clear VY
0141 %    fname = fullfile(rootpath, sub_ids{isub}, spmpath, 'SPM.mat');
0142 %    load(fname)
0143     
0144 %    filelist{isub} = strvcat(VY.fname);
0145     nvol = size(filelist{isub},1);    
0146     disp(sprintf('Loading data for subject: %s', sub_ids{isub}))
0147     for ivol = 1:nvol
0148       %
0149       %  Loop through all of the raw data files, grabbing the relevant voxels
0150       %
0151       vol = (spm_vol(deblank(filelist{isub}(ivol,:))));
0152       y(:,ivol,isub) = spm_sample_vol(vol, D(1,:),D(2,:),D(3,:),0);
0153     end
0154   end  % for isub=
0155 end