calculateDecodePerformance.m
2095645b
 % function [decodePerformance rawTimecourse ] = calculateDecodePerformance(des,timeLineStart, timeLineEnd, decodeDuration, svmargs, conditionList, sessionList, voxelList, classList, labelMap,normalize)
84376774
 function outputStruct = calculateDecodePerformance(inputStruct,SubjectID)
2095645b
 
 addpath 'libsvm-mat-2.88-1';
 
 outputStruct = struct;
 
84376774
 des             = inputStruct.(SubjectID).des;
2095645b
 timeLineStart   = inputStruct.frameShiftStart;
 timeLineEnd     = inputStruct.frameShiftEnd;
 decodeDuration  = inputStruct.decodeDuration;
 svmargs         = inputStruct.svmargs;
 sessionList     = inputStruct.sessionList;
84376774
 voxelList       = inputStruct.(SubjectID).voxelList;
2095645b
 % classList       = inputStruct.classList;
 % labelMap        = inputStruct.labelMap;
291ee33d
 smoothed       = inputStruct.smoothed;
2095645b
 globalStart     = inputStruct.psthStart;
 globalEnd       = inputStruct.psthEnd;
 baselineStart   = inputStruct.baselineStart;
 baselineEnd     = inputStruct.baselineEnd;
 eventList       = inputStruct.eventList;
 
 
 minPerformance = inf;
 maxPerformance = -inf;
 
 
         
         %Pro Voxel PSTH TIMELINE berechnen.
         %   timeshift mit pst-timeline durchf�hren.
         % psth-timeline -25 bis +15 zu RES Onset.
         
 %         eventList       = [9,11,13;10,12,14];
 %         globalStart     = -25;
 %         globalEnd       = 15;
 %         baselineStart   = -22;
 %         baselineEnd     = -20;
         
         
         for voxel = 1:size(voxelList,1)  % [[x;x],[y;y],[z;z]]
291ee33d
                 extr  = calculateImageData(voxelList(voxel,:),des,smoothed);
2095645b
                 rawdata=cell2mat({extr.mean}); % Raw Data
                 pst{voxel}  = calculatePST(des,globalStart,baselineStart,baselineEnd,globalEnd,eventList,rawdata,sessionList);
         end
 
         decodePerformance = [];
 
         for timeShift   = timeLineStart:1:timeLineEnd
             frameStart  = floor(-globalStart+1+timeShift - 0.5*decodeDuration);
             frameEnd    = min(ceil(frameStart+decodeDuration + 0.5*decodeDuration),-globalStart+globalEnd);
             
             tmp =[];
             anyvoxel = 1;
             for label = 1:size(pst{1,anyvoxel},2) 
                 for dp = 1:size(pst{1,anyvoxel}{1,label},1) % data point
                 row = label;
                     for voxel = 1:size(pst,2)
                         row = [row, pst{1,voxel}{1,label}(dp,frameStart:frameEnd)]; % label,value,value
                     end
                 tmp  = [tmp; row];
                 end
             end 
         
             svmdata      = tmp(:,2:size(tmp,2));
             svmlabel     = tmp(:,1);
             performance  = svmtrain(svmlabel, svmdata, svmargs);
 
             minPerformance = min(minPerformance,performance);
             maxPerformance = max(maxPerformance,performance);
 
             decodePerformance = [decodePerformance; performance];
         end
         
         outputStruct.decodePerformance  = decodePerformance;
         outputStruct.svmdata            = svmdata;
         outputStruct.svmlabel           = svmlabel;
         outputStruct.rawTimeCourse      = pst;
         outputStruct.minPerformance     = minPerformance;
         outputStruct.maxPerformance     = maxPerformance;
 
 % display(sprintf('Min CrossVal Accuracy: %g%% \t Max CrossVal Accuracy: %g%%',minPerformance,maxPerformance));
 end
 
 
291ee33d
 function extr = calculateImageData(voxelList,des,smoothed)
2095645b
 
 dtype='PSTH';
 
 switch dtype 
     case 'PSTH'
         V=des.xY.VY;
     case 'betas'
         V=des.Vbeta;
 end;
291ee33d
 
 if (~smoothed)
   for z=1:length(V) % Change Drive Letter!\
       % D:....SUBJECTID\session\swfanders...
       % D:....SUBJECTID\session\wfanders...
       tmp = findstr(filesep,V(z).fname);
       V(z).fname = strcat(V(z).fname(1:tmp(end)),V(z).fname(tmp(end)+2:end));
   end;
 end
2095645b
 
 % rad = 0; % one voxel
 % opt = 1; % xyz coordinates [mm]
 
 
 vox = voxelList;
 nRoi = size(vox,1);
 
 nImg = numel(V);
 
 for k=1:nImg
 	extr(k) = struct(...
         'val',   repmat(NaN, [1 nRoi]),...
 		'mean',  repmat(NaN, [1 nRoi]),...
 		'sum',   repmat(NaN, [1 nRoi]),...
 		'nvx',   repmat(NaN, [1 nRoi]),...
 		'posmm', repmat(NaN, [3 nRoi]),...
 		'posvx', repmat(NaN, [3 nRoi]));
 
     roicenter = round(inv(V(k).mat)*[vox, ones(nRoi,1)]');
 
 	for l = 1:nRoi
 
 %         if rad==0
             x = roicenter(1,l);
             y = roicenter(2,l);
             z = roicenter(3,l);
 %         else
 %             tmp = spm_imatrix(V(k).mat);
 %             vdim = tmp(7:9);
 %             vxrad = ceil((rad*ones(1,3))./(ones(nRoi,1)*vdim))';
 %             [x y z] = ndgrid(-vxrad(1,l):sign(vdim(1)):vxrad(1,l), ...
 %                       -vxrad(2,l):sign(vdim(2)):vxrad(2,l), ...
 %                       -vxrad(3,l):sign(vdim(3)):vxrad(3,l));
 %             sel = (x./vxrad(1,l)).^2 + (y./vxrad(2,l)).^2 + ...
 %                   (z./vxrad(3,l)).^2 <= 1;
 %             x = roicenter(1,l)+x(sel(:));
 %             y = roicenter(2,l)+y(sel(:));
 %             z = roicenter(3,l)+z(sel(:));
 %         end;
 		dat                 = spm_sample_vol(V(k), x, y, z,0);
 		[maxv maxi]         = max(dat);
 		tmp                 = V(k).mat*[x(maxi); y(maxi); z(maxi);1]; % Max Pos
 		extr(k).val(l)      = maxv;
 		extr(k).sum(l)      = sum(dat);
 		extr(k).mean(l)     = nanmean(dat);
         extr(k).nvx(l)      = numel(dat);
 		extr(k).posmm(:,l)  = tmp(1:3);
 		extr(k).posvx(:,l)  = [x(maxi); y(maxi); z(maxi)]; % Max Pos
 	end;
 
 end;
 end
 
 % disp(sprintf('Extracted at %.1f %.1f %.1f [xyz(mm)], average of %i voxel(s) [%.1fmm radius Sphere]',vox,length(x),rad));
 
 function pst = calculatePST(des,globalStart,baselineStart,baselineEnd,globalEnd,eventList,data,sessionList)
     bstart          = baselineStart;
     bend            = baselineEnd;
     edur            = 12;
     pre             =  globalStart;
     post            =  globalEnd;
     res             = 1;
 
     normz           = 'file';
     pm              = 0;
 
     lsess           = getNumberOfScans(des);
     nSessions       = getNumberOfSessions(des);
     tr              = 2;
 
     [evntrow evntcol]=size(eventList);
     
 
     hsec=str2num(des.xsDes.High_pass_Filter(8:end-3)); % Highpass filter [sec] from SPM.mat
 
     if strcmp(des.xBF.UNITS,'secs')
         unitsecs=1;
     end;
 
     nScansPerSession=getNumberOfScans(des);
     %stime=[0:tr:max(nScansPerSession)*tr+post-tr]; % Stimulus time for raw data plot
     stime=0:tr:max(nScansPerSession)*tr+round(post/tr)*tr-tr; % Stimulus time for raw data plot
 
 
 
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     % RUN
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 
     % Digital Highpass
     Rp=0.5;
     Rs=20;
     NO=1;
     Wp=1/((1/2/tr)/(1/hsec));
     [B, A] = ellip(NO,Rp,Rs,Wp,'high');
 
     sdata(1:max(nScansPerSession)+round(post/tr),1:nSessions)=nan; % Open Data Matrix
     for z=1:nSessions % Fill Data Matrix sessionwise
         sdata(1:nScansPerSession(z),z)=data(sum(nScansPerSession(1:z))-nScansPerSession(z)+1:sum(nScansPerSession(1:z)))';
     end;
 %         usdata=sdata; % Keep unfiltered data
 
     sdatamean=nanmean(nanmean(sdata(:,:)));
     for z=1:nSessions
 %             X(:,z)=[1:1:max(nScansPerSession)]'; % #Volume
         sdata(1:nScansPerSession(z),z)=filtfilt(B,A,sdata(1:nScansPerSession(z),z)); %Filter Data (Highpass)
     end;
     sdata=sdata+sdatamean;
 
 
     %%%%Parametric Modulation Modus%%%%
     if pm %Find Parameters for Event of Interest
         [imods modss mods erow evntrow eventList] = getParametricMappingEvents(eventList,evntrow,des,pmf);
     end;
     %%%%PM%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 
     for zr=1:evntrow
         n{zr}=0;
         nn{zr}=0; 
         nnn{zr}=0;
         sstart{zr}=1;
     end;
 
 
     sesst0=0; 
     for sessionID=sessionList
         if sessionID>1
             sesst0(sessionID)=sum(lsess(1:sessionID-1))*tr;  
         end;
         for zr=1:evntrow  %LABEL NUMBER, EVENT GROUP
             sstart{zr}=n{zr}+1;
             for ze=1:evntcol % EVENT INDEX in EventList
                 if ze==1 || (ze>1 && eventList(zr,ze)~=eventList(zr,ze-1))
                     for zz=1:length(des.Sess(sessionID).U(eventList(zr,ze)).ons) % EVENT REPETITION NUMBER
                         if ~unitsecs
                             des.Sess(sessionID).U(eventList(zr,ze)).ons(zz)=(des.Sess(sessionID).U(eventList(zr,ze)).ons(zz)-1)*tr;
                             des.Sess(sessionID).U(eventList(zr,ze)).dur(zz)=(des.Sess(sessionID).U(eventList(zr,ze)).dur(zz)-1)*tr;
                         end;
 
                         nnn{zr}=nnn{zr}+1; % INFO for rawdataplot start
                         if des.Sess(sessionID).U(eventList(zr,ze)).dur(zz)<edur
                             mev{zr}(nnn{zr},1:2)=[des.Sess(sessionID).U(eventList(zr,ze)).ons(zz)+sesst0(sessionID) edur]; % modeled event [onset length]
                         else
                             mev{zr}(nnn{zr},1:2)=[des.Sess(sessionID).U(eventList(zr,ze)).ons(zz)+sesst0(sessionID) des.Sess(sessionID).U(eventList(zr,ze)).dur(zz)];
                         end; % INFO for rawdataplot end
 
                         n{zr}=n{zr}+1;
                         pst{zr}(n{zr},:)=interp1(stime,sdata(:,sessionID),[des.Sess(sessionID).U(eventList(zr,ze)).ons(zz)+pre:res:des.Sess(sessionID).U(eventList(zr,ze)).ons(zz)+post],'linear');
                         if strcmp(normz,'epoc')
                             bline=nanmean(pst{zr}(n{zr},round(-pre/res+(bstart)/res+1):round(-pre/res+(bend)/res+1)));
                             if isnan(bline)
                                 pst{zr}(n{zr},1:-pre/res+post/res+1)=nan;
                             else
 %                                     nn{zr}=nn{zr}+1;
                                 pst{zr}(n{zr},:)=(pst{zr}(n{zr},:)-bline)/bline*100; % 'epoch-based' normalization
                             end;
                         end;
                     end;
                 end;
             end;
             if ~strcmp(normz,'epoc')
                 bline(zr)=nanmean(nanmean(pst{zr}(sstart{zr}:n{zr},-pre/res+(bstart)/res+1:-pre/res+(bend)/res+1)));
                 bstd(zr)=nanmean(nanstd(pst{zr}(sstart{zr}:n{zr},-pre/res+(bstart)/res+1:-pre/res+(bend)/res+1)));
                 nn{zr}=n{zr};
             end;
         end;
         if strcmp(normz,'filz')
             for zr=1:evntrow
                 pst{zr}(sstart{zr}:n{zr},:)=(pst{zr}(sstart{zr}:n{zr},:)-mean(bline))/mean(bstd); % session-based z-score normalization
             end;
         elseif strcmp(normz,'file')
             for zr=1:evntrow
                 pst{zr}(sstart{zr}:n{zr},:)=(pst{zr}(sstart{zr}:n{zr},:)-mean(bline))/mean(bline)*100; % session-based normalization
             end;
         end;
     end;
 end