function [seeds,rvalues] = analyzePRFcomputesupergridseeds(... res,stimulus,data,modelfun,maxpolydeg,dimdata,dimtime,typicalgain,noisereg,modeltype) % function [seeds,rvalues] = analyzePRF_computesupergridseeds(res,stimulus,data,modelfun,maxpolydeg,dimdata,dimtime,typicalgain,noisereg) % % is [R C] with the resolution of the stimuli % is a cell vector of time x (pixels+1) % is a cell vector of X x Y x Z x time (or XYZ x time) % is a function that accepts parameters (pp) and stimuli (dd) and outputs predicted time-series (time x 1) % is a vector of degrees (one element for each run) % is number of dimensions that pertain to voxels % is the dimension that is the time dimension % is a typical value for the gain in each time-series % is [] or a set of noise regressors (cell vector of matrices) % % this is an internal function called by analyzePRF.m. this function returns % as a matrix of dimensions X x Y x Z x parameters (or XYZ x parameters) % with the best seed from the super-grid. also, returns as X x Y x Z % (or XYZ x 1) with the corresponding correlation (r) values. % % history: % 2015/02/07 - make less memory intensive % internal notes: % - note that the gain seed is fake (it is not set the correct value but instead % to the ) % internal constants 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]; angs = linspacecircular(0,2*pi,16); expts = [0.5 0.25 0.125]; sdratio = [1 1.5 2 2.5 3]; normamp = [0 0.1 0.25 0.5 1]; % calc numvxs = prod(sizefull(data{1},dimdata)); % total number of voxels %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% calculate a long list of potential seeds % calc resmx = max(res); % calculate sigma gridding (pRF size is 1 px, 2 px, 4 px, ..., up to resmx) maxn = floor(log2(resmx)); ssindices = 2.^(0:maxn); % construct full list of seeds (seeds x params) [R C S G N] fprintf('constructing seeds.\n'); if strcmp(modeltype,'linear_hrf') || strcmp(modeltype,'linear_ephys') allseeds = zeros(length(eccs)*length(angs)*length(ssindices),4); elseif strcmp(modeltype,'css_hrf') || strcmp(modeltype,'css_ephys') allseeds = zeros(length(eccs)*length(angs)*length(ssindices)*length(expts),5); elseif strcmp(modeltype,'dog_hrf') || strcmp(modeltype,'dog_ephys') allseeds = zeros(length(eccs)*length(angs)*length(ssindices)*length(sdratio)*length(normamp),6); end cnt = 1; for p=1:length(eccs) for q=1:length(angs) if p==1 && q>1 % for the center-of-gaze, only do the first angle continue; end for s=1:length(ssindices) if strcmp(modeltype,'linear_hrf') || strcmp(modeltype,'linear_ephys') allseeds(cnt,:) = [(1+res(1))/2 - sin(angs(q)) * (eccs(p)*resmx) ... (1+res(2))/2 + cos(angs(q)) * (eccs(p)*resmx) ... ssindices(s) 1 ]; cnt = cnt + 1; elseif strcmp(modeltype,'css_hrf') || strcmp(modeltype,'css_ephys') for r=1:length(expts) allseeds(cnt,:) = [(1+res(1))/2 - sin(angs(q)) * (eccs(p)*resmx) ... (1+res(2))/2 + cos(angs(q)) * (eccs(p)*resmx) ... ssindices(s)*sqrt(expts(r)) 1 expts(r)]; cnt = cnt + 1; end elseif strcmp(modeltype,'dog_hrf') || strcmp(modeltype,'dog_ephys') for v=1:length(sdratio) for w=1:length(normamp) allseeds(cnt,:) = [(1+res(1))/2 - sin(angs(q)) * (eccs(p)*resmx) ... (1+res(2))/2 + cos(angs(q)) * (eccs(p)*resmx) ... ssindices(s) 1 sdratio(v) normamp(w)]; cnt = cnt + 1; end end end end end end allseeds(cnt:end,:) = []; % chop because of the omission above %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% generate the predicted time-series for each seed % generate predicted time-series [note that some time-series are all 0] predts = zeros(sum(cellfun(@(x) size(x,1),stimulus)),size(allseeds,1),'single'); % time x seeds temp = catcell(1,stimulus); fprintf('generating super-grid time-series...'); tic parfor p=1:size(allseeds,1) predts(:,p) = modelfun(allseeds(p,:),temp); end fprintf('done.'); toc clear temp; % % inspect for sanity on range and fineness [OBSOLETE] % figure; % for r=1:4:length(rrindices) % for c=1:4:length(ccindices) % for s=1:length(ssindices) % cnt = (r-1)*(length(ccindices)*length(ssindices)) + (c-1)*length(ssindices) + s; % clf; % plot(predts(:,cnt)'); % title(sprintf('r=%d, c=%d, s=%d',r,c,s)); % pause; % end % end % end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% prepare data and model predictions % construct polynomials, noise regressors, and projection matrix pregressors = {}; for p=1:length(maxpolydeg) pregressors{p} = constructpolynomialmatrix(size(data{p},dimtime),0:maxpolydeg(p)); if ~isempty(noisereg) pregressors{p} = cat(2,pregressors{p},noisereg{p}); end end pmatrix = projectionmatrix(blkdiag(pregressors{:})); % project out and scale to unit length predts = unitlength(pmatrix*predts, 1,[],0); % time x seeds [NOTE: some are all NaN] % OLD: datats = unitlength(pmatrix*squish(catcell(dimtime,data),dimdata)',1,[],0); % time x voxels %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% find the seed with the max correlation % compute correlation and find maximum for each voxel chunks = chunking(1:numvxs,100); rvalues = {}; bestseedix = {}; fprintf('finding best seed for each voxel.\n'); parfor p=1:length(chunks) % project out and scale to unit length datats = unitlength(pmatrix*catcell(2,cellfun(@(x) subscript(squish(x,dimdata),{chunks{p} ':'}),data,'UniformOutput',0))',1,[],0); % time x voxels % voxels x 1 with index of the best seed (max corr) [rvalues{p},bestseedix{p}] = max(datats'*predts,[],2); % voxels x seeds -> max corr along dim 2 [NaN is ok] end rvalues = catcell(1,rvalues); % voxels x 1 bestseedix = catcell(1,bestseedix); % voxels x 1 % prepare output rvalues = reshape(rvalues,[sizefull(data{1},dimdata) 1]); seeds = allseeds(bestseedix,:); % voxels x parameters seeds(:,4) = typicalgain; % set gain to typical gain seeds = reshape(seeds,[sizefull(data{1},dimdata) size(allseeds,2)]); % predts(:,p) = modelfun([allseeds(p,:) flatten(hrf)],temp); % is T x 1 with the HRF