% Example 1: Run analyzePRF on an example dataset
0001 %% Example 1: Run analyzePRF on an example dataset 0002 0003 %% Download dataset (if necessary) and add analyzePRF to the MATLAB path 0004 0005 setup; 0006 0007 %% Load in the data 0008 0009 % Load in the data 0010 load('exampledataset.mat'); 0011 0012 % Check the workspace 0013 whos 0014 %% 0015 0016 %% Inspect the data 0017 0018 % Check dimensionality of the data 0019 data 0020 %% 0021 0022 % There are four runs of data; each run consists of 150 time points (TR = 2 s). 0023 % The data have already been pre-processed for slice timing correction, motion 0024 % correction, and spatial undistortion. For simplicity, we have selected 0025 % 10 example voxels from the left hemisphere. Let's visualize the time-series 0026 % data for the second voxel. 0027 temp = cellfun(@(x) x(2,:),data,'UniformOutput',0); 0028 figure; hold on; 0029 set(gcf,'Units','points','Position',[100 100 600 150]); 0030 plot(cat(2,temp{:}),'r-'); 0031 straightline(150*(1:4)+.5,'v','g-'); 0032 xlabel('TR'); 0033 ylabel('BOLD signal'); 0034 ax = axis; 0035 axis([.5 600+.5 ax(3:4)]); 0036 title('Time-series data'); 0037 %% 0038 0039 %% Inspect the stimuli 0040 0041 % Check dimensionality of the stimuli 0042 stimulus 0043 %% 0044 0045 % The stimulus images have been prepared at a resolution of 100 pixels x 100 pixels. 0046 % There are 300 images per run because we have prepared the images at a time resolution 0047 % of 1 second. (Note that this is faster than the data sampling rate. When analyzing 0048 % the data, we will resample the data to match the stimulus rate.) Let's inspect a 0049 % few of the stimulus images in the first run. 0050 figure; 0051 set(gcf,'Units','points','Position',[100 100 700 300]); 0052 for p=1:3 0053 subplot(1,3,p); hold on; 0054 num = 239+2*p; 0055 imagesc(stimulus{1}(:,:,num),[0 1]); 0056 axis image tight; 0057 set(gca,'YDir','reverse'); 0058 colormap(gray); 0059 title(sprintf('Image number %d',num)); 0060 end 0061 %% 0062 0063 % Notice that the range of values is 0 to 1 (0 indicates that the gray background was 0064 % present; 1 indicates that the stimulus was present). For these stimulus images, 0065 % the stimulus is a bar that moves downward and to the left. 0066 0067 %% Analyze the data 0068 0069 % Start parallel MATLAB to speed up execution. 0070 if matlabpool('size')==0 0071 matlabpool open; 0072 end 0073 0074 % We need to resample the data to match the temporal rate of the stimulus. Here we use 0075 % cubic interpolation to increase the rate of the data from 2 seconds to 1 second (note 0076 % that the first time point is preserved and the total data duration stays the same). 0077 data = tseriesinterp(data,2,1,2); 0078 0079 % Finally, we analyze the data using analyzePRF. The third argument is the TR, which 0080 % is now 1 second. The default analysis strategy involves two generic initial seeds 0081 % and an initial seed that is computed by performing a grid search. This last seed is 0082 % a little costly to compute, so to save time, we set 'seedmode' to [0 1] which means 0083 % to just use the two generic initial seeds. We suppress command-window output by 0084 % setting 'display' to 'off'. 0085 results = analyzePRF(stimulus,data,1,struct('seedmode',[0 1],'display','off')); 0086 %% 0087 0088 % Note that because of the use of parfor, the command-window output for different 0089 % voxels will come in at different times (and so the text above is not really 0090 % readable). 0091 0092 %% Inspect the results 0093 0094 % The stimulus is 100 pixels (in both height and weight), and this corresponds to 0095 % 10 degrees of visual angle. To convert from pixels to degreees, we multiply 0096 % by 10/100. 0097 cfactor = 10/100; 0098 0099 % Visualize the location of each voxel's pRF 0100 figure; hold on; 0101 set(gcf,'Units','points','Position',[100 100 400 400]); 0102 cmap = jet(size(results.ang,1)); 0103 for p=1:size(results.ang,1) 0104 xpos = results.ecc(p) * cos(results.ang(p)/180*pi) * cfactor; 0105 ypos = results.ecc(p) * sin(results.ang(p)/180*pi) * cfactor; 0106 ang = results.ang(p)/180*pi; 0107 sd = results.rfsize(p) * cfactor; 0108 h = drawellipse(xpos,ypos,ang,2*sd,2*sd); % circle at +/- 2 pRF sizes 0109 set(h,'Color',cmap(p,:),'LineWidth',2); 0110 set(scatter(xpos,ypos,'r.'),'CData',cmap(p,:)); 0111 end 0112 drawrectangle(0,0,10,10,'k-'); % square indicating stimulus extent 0113 axis([-10 10 -10 10]); 0114 straightline(0,'h','k-'); % line indicating horizontal meridian 0115 straightline(0,'v','k-'); % line indicating vertical meridian 0116 axis square; 0117 set(gca,'XTick',-10:2:10,'YTick',-10:2:10); 0118 xlabel('X-position (deg)'); 0119 ylabel('Y-position (deg)'); 0120 %% 0121 0122 % Please see the example2.m script for an example of how to inspect the model fit 0123 % and compare it to the data.