function seeds = analyzePRF_computesupergridseeds(res,stimulus,data,modelfun,maxpolydeg,dimdata,dimtime,typicalgain) <res> is [R C] with the resolution of the stimuli <stimulus> is a cell vector of time x (pixels+1) <data> is a cell vector of X x Y x Z x time (or XYZ x time) <modelfun> is a function that accepts parameters (pp) and stimuli (dd) and outputs predicted time-series (time x 1) <maxpolydeg> is a vector of degrees (one element for each run) <dimdata> is number of dimensions that pertain to voxels <dimtime> is the dimension that is the time dimension <typicalgain> is a typical value for the gain in each time-series this is an internal function called by analyzePRF.m. this function returns <seeds> as a matrix of dimensions X x Y x Z x parameters (or XYZ x parameters) with the best seed from the super-grid.
0001 function seeds = analyzePRF_computesupergridseeds(res,stimulus,data,modelfun,maxpolydeg,dimdata,dimtime,typicalgain) 0002 0003 % function seeds = analyzePRF_computesupergridseeds(res,stimulus,data,modelfun,maxpolydeg,dimdata,dimtime,typicalgain) 0004 % 0005 % <res> is [R C] with the resolution of the stimuli 0006 % <stimulus> is a cell vector of time x (pixels+1) 0007 % <data> is a cell vector of X x Y x Z x time (or XYZ x time) 0008 % <modelfun> is a function that accepts parameters (pp) and stimuli (dd) and outputs predicted time-series (time x 1) 0009 % <maxpolydeg> is a vector of degrees (one element for each run) 0010 % <dimdata> is number of dimensions that pertain to voxels 0011 % <dimtime> is the dimension that is the time dimension 0012 % <typicalgain> is a typical value for the gain in each time-series 0013 % 0014 % this is an internal function called by analyzePRF.m. this function returns <seeds> 0015 % as a matrix of dimensions X x Y x Z x parameters (or XYZ x parameters) 0016 % with the best seed from the super-grid. 0017 0018 % internal notes: 0019 % - note that the gain seed is fake (it is not set the correct value but instead 0020 % to the <typicalgain>) 0021 0022 % internal constants 0023 eccs = [0 0.00551 0.014 0.0269 0.0459 0.0731 0.112 0.166 0.242 0.348 0.498 0.707 1]; 0024 angs = linspacecircular(0,2*pi,16); 0025 expts = [0.5 0.25 0.125]; 0026 0027 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% calculate a long list of potential seeds 0028 0029 % calc 0030 resmx = max(res); 0031 0032 % calculate sigma gridding (pRF size is 1 px, 2 px, 4 px, ..., up to resmx) 0033 maxn = floor(log2(resmx)); 0034 ssindices = 2.^(0:maxn); 0035 0036 % construct full list of seeds (seeds x params) [R C S G N] 0037 fprintf('constructing seeds.\n'); 0038 allseeds = zeros(length(eccs)*length(angs)*length(ssindices)*length(expts),5); 0039 cnt = 1; 0040 for p=1:length(eccs) 0041 for q=1:length(angs) 0042 if p==1 && q>1 % for the center-of-gaze, only do the first angle 0043 continue; 0044 end 0045 for s=1:length(ssindices) 0046 for r=1:length(expts) 0047 allseeds(cnt,:) = [(1+res(1))/2 - sin(angs(q)) * (eccs(p)*resmx) ... 0048 (1+res(2))/2 + cos(angs(q)) * (eccs(p)*resmx) ... 0049 ssindices(s)*sqrt(expts(r)) 1 expts(r)]; 0050 cnt = cnt + 1; 0051 end 0052 end 0053 end 0054 end 0055 allseeds(cnt:end,:) = []; % chop because of the omission above 0056 0057 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% generate the predicted time-series for each seed 0058 0059 % generate predicted time-series [note that some time-series are all 0] 0060 predts = zeros(sum(cellfun(@(x) size(x,1),stimulus)),size(allseeds,1),'single'); % time x seeds 0061 temp = catcell(1,stimulus); 0062 fprintf('generating super-grid time-series...'); tic 0063 parfor p=1:size(allseeds,1) 0064 predts(:,p) = modelfun(allseeds(p,:),temp); 0065 end 0066 fprintf('done.'); toc 0067 clear temp; 0068 0069 % % inspect for sanity on range and fineness [OBSOLETE] 0070 % figure; 0071 % for r=1:4:length(rrindices) 0072 % for c=1:4:length(ccindices) 0073 % for s=1:length(ssindices) 0074 % cnt = (r-1)*(length(ccindices)*length(ssindices)) + (c-1)*length(ssindices) + s; 0075 % clf; 0076 % plot(predts(:,cnt)'); 0077 % title(sprintf('r=%d, c=%d, s=%d',r,c,s)); 0078 % pause; 0079 % end 0080 % end 0081 % end 0082 0083 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% prepare data and model predictions 0084 0085 % construct polynomials and projection matrix 0086 polyregressors = {}; 0087 for p=1:length(maxpolydeg) 0088 polyregressors{p} = constructpolynomialmatrix(size(data{p},dimtime),0:maxpolydeg(p)); 0089 end 0090 pmatrix = projectionmatrix(blkdiag(polyregressors{:})); 0091 0092 % project out polynomials and scale to unit length 0093 predts = unitlength(pmatrix*predts, 1,[],0); % time x seeds [NOTE: some are all NaN] 0094 datats = unitlength(pmatrix*squish(catcell(dimtime,data),dimdata)',1,[],0); % time x voxels 0095 0096 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% find the seed with the max correlation 0097 0098 % compute correlation and find maximum for each voxel 0099 chunks = chunking(1:size(datats,2),100); 0100 bestseedix = {}; 0101 fprintf('finding best seed for each voxel.\n'); 0102 parfor p=1:length(chunks) 0103 % voxels x 1 with index of the best seed (max corr) 0104 [mx,bestseedix{p}] = max(datats(:,chunks{p})' * predts,[],2); % voxels x seeds -> max corr along dim 2 [NaN is ok] 0105 end 0106 bestseedix = catcell(1,bestseedix); % voxels x 1 0107 0108 % prepare output 0109 seeds = allseeds(bestseedix,:); % voxels x parameters 0110 seeds(:,4) = typicalgain; % set gain to typical gain 0111 seeds = reshape(seeds,[sizefull(data{1},dimdata) size(allseeds,2)]); 0112 0113 0114 0115 0116 0117 % predts(:,p) = modelfun([allseeds(p,:) flatten(hrf)],temp); 0118 % <hrf> is T x 1 with the HRF