Browse code

moving things to private, some modifications

git-svn-id: https://svn.discofish.de/MATLAB/spmtoolbox/SVMCrossVal@110 83ab2cfd-5345-466c-8aeb-2b2739fb922d

Christoph Budziszewski authored on16/01/2009 02:48:58
Showing1 changed files
1 1
deleted file mode 100644
... ...
@@ -1,95 +0,0 @@
1
-% function [decodePerformance rawTimecourse ] = calculateDecodePerformance(des,timeLineStart, timeLineEnd, decodeDuration, svmargs, conditionList, sessionList, voxelList, classList, labelMap,normalize)
2
-function outputStruct = calculateDecodePerformance(inputStruct,SubjectID)
3
-
4
-addpath 'libsvm-mat-2.88-1';
5
-
6
-SINGLE = 1;
7
-
8
-outputStruct = struct;
9
-
10
-namehelper      = strcat('s',SubjectID);
11
-des             = inputStruct.(namehelper).des;
12
-timeLineStart   = inputStruct.frameShiftStart;
13
-timeLineEnd     = inputStruct.frameShiftEnd;
14
-decodeDuration  = inputStruct.decodeDuration;
15
-svmargs         = inputStruct.svmargs;
16
-sessionList     = inputStruct.sessionList;
17
-voxelList       = inputStruct.(namehelper).voxelList;
18
-% classList       = inputStruct.classList;
19
-% labelMap        = inputStruct.labelMap;
20
-smoothed       = inputStruct.smoothed;
21
-globalStart     = inputStruct.psthStart;
22
-globalEnd       = inputStruct.psthEnd;
23
-baselineStart   = inputStruct.baselineStart;
24
-baselineEnd     = inputStruct.baselineEnd;
25
-eventList       = inputStruct.eventList;
26
-labelMap        = inputStruct.labelMap;
27
-
28
-
29
-minPerformance = inf;
30
-maxPerformance = -inf;
31
-
32
-
33
-        
34
-        %Pro Voxel PSTH TIMELINE berechnen.
35
-        %   timeshift mit pst-timeline durchf�hren.
36
-        % psth-timeline -25 bis +15 zu RES Onset.
37
-        
38
-%         eventList       = [9,11,13;10,12,14];
39
-%         globalStart     = -25;
40
-%         globalEnd       = 15;
41
-%         baselineStart   = -22;
42
-%         baselineEnd     = -20;
43
-        
44
-        
45
-        for voxel = 1:size(voxelList,1)  % [[x;x],[y;y],[z;z]]
46
-                extr  = calculateImageData(voxelList(voxel,:),des,smoothed);
47
-                rawdata=cell2mat({extr.mean}); % Raw Data
48
-                pst{voxel}  = calculatePST(des,globalStart,baselineStart,baselineEnd,globalEnd,eventList,rawdata,sessionList);
49
-        end
50
-        
51
-        decodePerformance = [];
52
-
53
-        for timeShift   = timeLineStart:1:timeLineEnd
54
-            frameStart  = floor(-globalStart+1+timeShift - 0.5*decodeDuration);
55
-            frameEnd    = min(ceil(frameStart+decodeDuration + 0.5*decodeDuration),-globalStart+globalEnd);
56
-            
57
-            tmp =[];
58
-            anyvoxel = 1;
59
-            for pstConditionGroup = 1:size(pst{1,anyvoxel},2) 
60
-                for dp = 1:size(pst{1,anyvoxel}{1,pstConditionGroup},1) % data point
61
-                  row = getSVMLabel(labelMap,eventList(pstConditionGroup,1));
62
-                    for voxel = 1:size(pst,2)
63
-                        row = [row, pst{1,voxel}{1,pstConditionGroup}(dp,frameStart:frameEnd)]; % label,value,value
64
-                    end
65
-                tmp  = [tmp; row];
66
-                end
67
-            end 
68
-        
69
-            svmdata      = tmp(:,2:size(tmp,2));
70
-            svmlabel     = tmp(:,1);
71
-            
72
-%             RANDOMIZE INPUT
73
-%             rndindex  = randperm(length(svmlabel));
74
-%             svmdata   = svmdata(rndindex,:);
75
-%             svmlabel  = svmlabel(rndindex);
76
-            
77
-            if SINGLE
78
-                performance  = svmtrain(svmlabel, svmdata, svmargs);
79
-
80
-                minPerformance = min(minPerformance,performance);
81
-                maxPerformance = max(maxPerformance,performance);
82
-
83
-                decodePerformance = [decodePerformance; performance];
84
-            end
85
-            
86
-        end
87
-        
88
-        outputStruct.decodePerformance  = decodePerformance;
89
-        outputStruct.svmdata            = svmdata;
90
-        outputStruct.svmlabel           = svmlabel;
91
-        outputStruct.rawTimeCourse      = pst;
92
-        outputStruct.minPerformance     = minPerformance;
93
-        outputStruct.maxPerformance     = maxPerformance;
94
-end
95
-
Browse code

New Normalization Methods in CalculatePST

git-svn-id: https://svn.discofish.de/MATLAB/spmtoolbox/SVMCrossVal@104 83ab2cfd-5345-466c-8aeb-2b2739fb922d

Axel Lindner authored on09/01/2009 16:47:12
Showing1 changed files
... ...
@@ -69,6 +69,7 @@ maxPerformance = -inf;
69 69
             svmdata      = tmp(:,2:size(tmp,2));
70 70
             svmlabel     = tmp(:,1);
71 71
             
72
+%             RANDOMIZE INPUT
72 73
 %             rndindex  = randperm(length(svmlabel));
73 74
 %             svmdata   = svmdata(rndindex,:);
74 75
 %             svmlabel  = svmlabel(rndindex);
Browse code

plot optimiert. normalisierung pst, zwischenversion UI Defaults verbessert

git-svn-id: https://svn.discofish.de/MATLAB/spmtoolbox/SVMCrossVal@103 83ab2cfd-5345-466c-8aeb-2b2739fb922d

Axel Lindner authored on09/01/2009 15:52:07
Showing1 changed files
... ...
@@ -47,7 +47,7 @@ maxPerformance = -inf;
47 47
                 rawdata=cell2mat({extr.mean}); % Raw Data
48 48
                 pst{voxel}  = calculatePST(des,globalStart,baselineStart,baselineEnd,globalEnd,eventList,rawdata,sessionList);
49 49
         end
50
-
50
+        
51 51
         decodePerformance = [];
52 52
 
53 53
         for timeShift   = timeLineStart:1:timeLineEnd
... ...
@@ -58,7 +58,7 @@ maxPerformance = -inf;
58 58
             anyvoxel = 1;
59 59
             for pstConditionGroup = 1:size(pst{1,anyvoxel},2) 
60 60
                 for dp = 1:size(pst{1,anyvoxel}{1,pstConditionGroup},1) % data point
61
-                  row = getSVMLabel(labelMap,eventList(pstConditionGroup,1)) 
61
+                  row = getSVMLabel(labelMap,eventList(pstConditionGroup,1));
62 62
                     for voxel = 1:size(pst,2)
63 63
                         row = [row, pst{1,voxel}{1,pstConditionGroup}(dp,frameStart:frameEnd)]; % label,value,value
64 64
                     end
... ...
@@ -69,6 +69,10 @@ maxPerformance = -inf;
69 69
             svmdata      = tmp(:,2:size(tmp,2));
70 70
             svmlabel     = tmp(:,1);
71 71
             
72
+%             rndindex  = randperm(length(svmlabel));
73
+%             svmdata   = svmdata(rndindex,:);
74
+%             svmlabel  = svmlabel(rndindex);
75
+            
72 76
             if SINGLE
73 77
                 performance  = svmtrain(svmlabel, svmdata, svmargs);
74 78
 
Browse code

kleine änderungen code move

git-svn-id: https://svn.discofish.de/MATLAB/spmtoolbox/SVMCrossVal@102 83ab2cfd-5345-466c-8aeb-2b2739fb922d

Christoph Budziszewski authored on08/01/2009 11:44:25
Showing1 changed files
... ...
@@ -3,6 +3,8 @@ function outputStruct = calculateDecodePerformance(inputStruct,SubjectID)
3 3
 
4 4
 addpath 'libsvm-mat-2.88-1';
5 5
 
6
+SINGLE = 1;
7
+
6 8
 outputStruct = struct;
7 9
 
8 10
 namehelper      = strcat('s',SubjectID);
... ...
@@ -66,12 +68,16 @@ maxPerformance = -inf;
66 68
         
67 69
             svmdata      = tmp(:,2:size(tmp,2));
68 70
             svmlabel     = tmp(:,1);
69
-            performance  = svmtrain(svmlabel, svmdata, svmargs);
71
+            
72
+            if SINGLE
73
+                performance  = svmtrain(svmlabel, svmdata, svmargs);
70 74
 
71
-            minPerformance = min(minPerformance,performance);
72
-            maxPerformance = max(maxPerformance,performance);
75
+                minPerformance = min(minPerformance,performance);
76
+                maxPerformance = max(maxPerformance,performance);
73 77
 
74
-            decodePerformance = [decodePerformance; performance];
78
+                decodePerformance = [decodePerformance; performance];
79
+            end
80
+            
75 81
         end
76 82
         
77 83
         outputStruct.decodePerformance  = decodePerformance;
... ...
@@ -80,215 +86,5 @@ maxPerformance = -inf;
80 86
         outputStruct.rawTimeCourse      = pst;
81 87
         outputStruct.minPerformance     = minPerformance;
82 88
         outputStruct.maxPerformance     = maxPerformance;
83
-
84
-% display(sprintf('Min CrossVal Accuracy: %g%% \t Max CrossVal Accuracy: %g%%',minPerformance,maxPerformance));
85
-end
86
-
87
-
88
-function extr = calculateImageData(voxelList,des,smoothed)
89
-
90
-dtype='PSTH';
91
-
92
-switch dtype 
93
-    case 'PSTH'
94
-        V=des.xY.VY;
95
-    case 'betas'
96
-        V=des.Vbeta;
97
-end;
98
-
99
-if V(1).fname(1)~='D'
100
-     for z=1:length(V) % Change Drive Letter!\
101
-      	V(z).fname(1) = 'D';
102
-     end; 
103
-end
104
-
105
-if (~smoothed)
106
-  for z=1:length(V) % Change Drive Letter!\
107
-      % D:....SUBJECTID\session\swfanders...
108
-      % D:....SUBJECTID\session\wfanders...
109
-      tmp = findstr(filesep,V(z).fname);
110
-      V(z).fname = strcat(V(z).fname(1:tmp(end)),V(z).fname(tmp(end)+2:end));
111
-  end;
112
-end
113
-
114
-% rad = 0; % one voxel
115
-% opt = 1; % xyz coordinates [mm]
116
-
117
-
118
-vox = voxelList;
119
-nRoi = size(vox,1);
120
-
121
-nImg = numel(V);
122
-
123
-for k=1:nImg
124
-	extr(k) = struct(...
125
-        'val',   repmat(NaN, [1 nRoi]),...
126
-		'mean',  repmat(NaN, [1 nRoi]),...
127
-		'sum',   repmat(NaN, [1 nRoi]),...
128
-		'nvx',   repmat(NaN, [1 nRoi]),...
129
-		'posmm', repmat(NaN, [3 nRoi]),...
130
-		'posvx', repmat(NaN, [3 nRoi]));
131
-
132
-    roicenter = round(inv(V(k).mat)*[vox, ones(nRoi,1)]');
133
-
134
-	for l = 1:nRoi
135
-
136
-%         if rad==0
137
-            x = roicenter(1,l);
138
-            y = roicenter(2,l);
139
-            z = roicenter(3,l);
140
-%         else
141
-%             tmp = spm_imatrix(V(k).mat);
142
-%             vdim = tmp(7:9);
143
-%             vxrad = ceil((rad*ones(1,3))./(ones(nRoi,1)*vdim))';
144
-%             [x y z] = ndgrid(-vxrad(1,l):sign(vdim(1)):vxrad(1,l), ...
145
-%                       -vxrad(2,l):sign(vdim(2)):vxrad(2,l), ...
146
-%                       -vxrad(3,l):sign(vdim(3)):vxrad(3,l));
147
-%             sel = (x./vxrad(1,l)).^2 + (y./vxrad(2,l)).^2 + ...
148
-%                   (z./vxrad(3,l)).^2 <= 1;
149
-%             x = roicenter(1,l)+x(sel(:));
150
-%             y = roicenter(2,l)+y(sel(:));
151
-%             z = roicenter(3,l)+z(sel(:));
152
-%         end;
153
-
154
-
155
-		dat                 = spm_sample_vol(V(k), x, y, z,0);
156
-		[maxv maxi]         = max(dat);
157
-		tmp                 = V(k).mat*[x(maxi); y(maxi); z(maxi);1]; % Max Pos
158
-		extr(k).val(l)      = maxv;
159
-		extr(k).sum(l)      = sum(dat);
160
-		extr(k).mean(l)     = nanmean(dat);
161
-        extr(k).nvx(l)      = numel(dat);
162
-		extr(k).posmm(:,l)  = tmp(1:3);
163
-		extr(k).posvx(:,l)  = [x(maxi); y(maxi); z(maxi)]; % Max Pos
164
-	end;
165
-
166
-end;
167 89
 end
168 90
 
169
-% disp(sprintf('Extracted at %.1f %.1f %.1f [xyz(mm)], average of %i voxel(s) [%.1fmm radius Sphere]',vox,length(x),rad));
170
-
171
-function pst = calculatePST(des,globalStart,baselineStart,baselineEnd,globalEnd,eventList,data,sessionList)
172
-    bstart          = baselineStart;
173
-    bend            = baselineEnd;
174
-    edur            = 12;
175
-    pre             =  globalStart;
176
-    post            =  globalEnd;
177
-    res             = 1;
178
-
179
-    normz           = 'file';
180
-    pm              = 0;
181
-
182
-    lsess           = getNumberOfScans(des);
183
-    nSessions       = getNumberOfSessions(des);
184
-    tr              = 2;
185
-
186
-    [evntrow evntcol]=size(eventList);
187
-    
188
-
189
-    hsec=str2num(des.xsDes.High_pass_Filter(8:end-3)); % Highpass filter [sec] from SPM.mat
190
-
191
-    if strcmp(des.xBF.UNITS,'secs')
192
-        unitsecs=1;
193
-    end;
194
-
195
-    nScansPerSession=getNumberOfScans(des);
196
-    %stime=[0:tr:max(nScansPerSession)*tr+post-tr]; % Stimulus time for raw data plot
197
-    stime=0:tr:max(nScansPerSession)*tr+round(post/tr)*tr-tr; % Stimulus time for raw data plot
198
-
199
-
200
-
201
-    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
202
-    % RUN
203
-    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
204
-
205
-
206
-    % Digital Highpass
207
-    Rp=0.5;
208
-    Rs=20;
209
-    NO=1;
210
-    Wp=1/((1/2/tr)/(1/hsec));
211
-    [B, A] = ellip(NO,Rp,Rs,Wp,'high');
212
-
213
-    sdata(1:max(nScansPerSession)+round(post/tr),1:nSessions)=nan; % Open Data Matrix
214
-    for z=1:nSessions % Fill Data Matrix sessionwise
215
-        sdata(1:nScansPerSession(z),z)=data(sum(nScansPerSession(1:z))-nScansPerSession(z)+1:sum(nScansPerSession(1:z)))';
216
-    end;
217
-%         usdata=sdata; % Keep unfiltered data
218
-
219
-    sdatamean=nanmean(nanmean(sdata(:,:)));
220
-    for z=1:nSessions
221
-%             X(:,z)=[1:1:max(nScansPerSession)]'; % #Volume
222
-        sdata(1:nScansPerSession(z),z)=filtfilt(B,A,sdata(1:nScansPerSession(z),z)); %Filter Data (Highpass)
223
-    end;
224
-    sdata=sdata+sdatamean;
225
-
226
-
227
-    %%%%Parametric Modulation Modus%%%%
228
-    if pm %Find Parameters for Event of Interest
229
-        [imods modss mods erow evntrow eventList] = getParametricMappingEvents(eventList,evntrow,des,pmf);
230
-    end;
231
-    %%%%PM%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
232
-
233
-
234
-    for zr=1:evntrow
235
-        n{zr}=0;
236
-        nn{zr}=0; 
237
-        nnn{zr}=0;
238
-        sstart{zr}=1;
239
-    end;
240
-
241
-
242
-    sesst0=0; 
243
-    for sessionID=sessionList
244
-        if sessionID>1
245
-            sesst0(sessionID)=sum(lsess(1:sessionID-1))*tr;  
246
-        end;
247
-        for zr=1:evntrow  %LABEL NUMBER, EVENT GROUP
248
-            sstart{zr}=n{zr}+1;
249
-            for ze=1:evntcol % EVENT INDEX in EventList
250
-                if ze==1 || (ze>1 && eventList(zr,ze)~=eventList(zr,ze-1))
251
-                    for zz=1:length(des.Sess(sessionID).U(eventList(zr,ze)).ons) % EVENT REPETITION NUMBER
252
-                        if ~unitsecs
253
-                            des.Sess(sessionID).U(eventList(zr,ze)).ons(zz)=(des.Sess(sessionID).U(eventList(zr,ze)).ons(zz)-1)*tr;
254
-                            des.Sess(sessionID).U(eventList(zr,ze)).dur(zz)=(des.Sess(sessionID).U(eventList(zr,ze)).dur(zz)-1)*tr;
255
-                        end;
256
-
257
-                        nnn{zr}=nnn{zr}+1; % INFO for rawdataplot start
258
-                        if des.Sess(sessionID).U(eventList(zr,ze)).dur(zz)<edur
259
-                            mev{zr}(nnn{zr},1:2)=[des.Sess(sessionID).U(eventList(zr,ze)).ons(zz)+sesst0(sessionID) edur]; % modeled event [onset length]
260
-                        else
261
-                            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)];
262
-                        end; % INFO for rawdataplot end
263
-
264
-                        n{zr}=n{zr}+1;
265
-                        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');
266
-                        if strcmp(normz,'epoc')
267
-                            bline=nanmean(pst{zr}(n{zr},round(-pre/res+(bstart)/res+1):round(-pre/res+(bend)/res+1)));
268
-                            if isnan(bline)
269
-                                pst{zr}(n{zr},1:-pre/res+post/res+1)=nan;
270
-                            else
271
-%                                     nn{zr}=nn{zr}+1;
272
-                                pst{zr}(n{zr},:)=(pst{zr}(n{zr},:)-bline)/bline*100; % 'epoch-based' normalization
273
-                            end;
274
-                        end;
275
-                    end;
276
-                end;
277
-            end;
278
-            if ~strcmp(normz,'epoc')
279
-                bline(zr)=nanmean(nanmean(pst{zr}(sstart{zr}:n{zr},-pre/res+(bstart)/res+1:-pre/res+(bend)/res+1)));
280
-                bstd(zr)=nanmean(nanstd(pst{zr}(sstart{zr}:n{zr},-pre/res+(bstart)/res+1:-pre/res+(bend)/res+1)));
281
-                nn{zr}=n{zr};
282
-            end;
283
-        end;
284
-        if strcmp(normz,'filz')
285
-            for zr=1:evntrow
286
-                pst{zr}(sstart{zr}:n{zr},:)=(pst{zr}(sstart{zr}:n{zr},:)-mean(bline))/mean(bstd); % session-based z-score normalization
287
-            end;
288
-        elseif strcmp(normz,'file')
289
-            for zr=1:evntrow
290
-                pst{zr}(sstart{zr}:n{zr},:)=(pst{zr}(sstart{zr}:n{zr},:)-mean(bline))/mean(bline)*100; % session-based normalization
291
-            end;
292
-        end;
293
-    end;
294
-end
Browse code

new LabelMap new svm grouping method

git-svn-id: https://svn.discofish.de/MATLAB/spmtoolbox/SVMCrossVal@100 83ab2cfd-5345-466c-8aeb-2b2739fb922d

Christoph Budziszewski authored on07/01/2009 18:26:33
Showing1 changed files
... ...
@@ -21,6 +21,7 @@ globalEnd       = inputStruct.psthEnd;
21 21
 baselineStart   = inputStruct.baselineStart;
22 22
 baselineEnd     = inputStruct.baselineEnd;
23 23
 eventList       = inputStruct.eventList;
24
+labelMap        = inputStruct.labelMap;
24 25
 
25 26
 
26 27
 minPerformance = inf;
... ...
@@ -53,11 +54,11 @@ maxPerformance = -inf;
53 54
             
54 55
             tmp =[];
55 56
             anyvoxel = 1;
56
-            for label = 1:size(pst{1,anyvoxel},2) 
57
-                for dp = 1:size(pst{1,anyvoxel}{1,label},1) % data point
58
-                row = label;
57
+            for pstConditionGroup = 1:size(pst{1,anyvoxel},2) 
58
+                for dp = 1:size(pst{1,anyvoxel}{1,pstConditionGroup},1) % data point
59
+                  row = getSVMLabel(labelMap,eventList(pstConditionGroup,1)) 
59 60
                     for voxel = 1:size(pst,2)
60
-                        row = [row, pst{1,voxel}{1,label}(dp,frameStart:frameEnd)]; % label,value,value
61
+                        row = [row, pst{1,voxel}{1,pstConditionGroup}(dp,frameStart:frameEnd)]; % label,value,value
61 62
                     end
62 63
                 tmp  = [tmp; row];
63 64
                 end
Browse code

New5 Studie, Var-Name ersetzungs bug

git-svn-id: https://svn.discofish.de/MATLAB/spmtoolbox/SVMCrossVal@98 83ab2cfd-5345-466c-8aeb-2b2739fb922d

Christoph Budziszewski authored on07/01/2009 12:57:14
Showing1 changed files
... ...
@@ -5,13 +5,14 @@ addpath 'libsvm-mat-2.88-1';
5 5
 
6 6
 outputStruct = struct;
7 7
 
8
-des             = inputStruct.(SubjectID).des;
8
+namehelper      = strcat('s',SubjectID);
9
+des             = inputStruct.(namehelper).des;
9 10
 timeLineStart   = inputStruct.frameShiftStart;
10 11
 timeLineEnd     = inputStruct.frameShiftEnd;
11 12
 decodeDuration  = inputStruct.decodeDuration;
12 13
 svmargs         = inputStruct.svmargs;
13 14
 sessionList     = inputStruct.sessionList;
14
-voxelList       = inputStruct.(SubjectID).voxelList;
15
+voxelList       = inputStruct.(namehelper).voxelList;
15 16
 % classList       = inputStruct.classList;
16 17
 % labelMap        = inputStruct.labelMap;
17 18
 smoothed       = inputStruct.smoothed;
... ...
@@ -94,6 +95,12 @@ switch dtype
94 95
         V=des.Vbeta;
95 96
 end;
96 97
 
98
+if V(1).fname(1)~='D'
99
+     for z=1:length(V) % Change Drive Letter!\
100
+      	V(z).fname(1) = 'D';
101
+     end; 
102
+end
103
+
97 104
 if (~smoothed)
98 105
   for z=1:length(V) % Change Drive Letter!\
99 106
       % D:....SUBJECTID\session\swfanders...
... ...
@@ -142,6 +149,8 @@ for k=1:nImg
142 149
 %             y = roicenter(2,l)+y(sel(:));
143 150
 %             z = roicenter(3,l)+z(sel(:));
144 151
 %         end;
152
+
153
+
145 154
 		dat                 = spm_sample_vol(V(k), x, y, z,0);
146 155
 		[maxv maxi]         = max(dat);
147 156
 		tmp                 = V(k).mat*[x(maxi); y(maxi); z(maxi);1]; % Max Pos
Browse code

multi-subject support

git-svn-id: https://svn.discofish.de/MATLAB/spmtoolbox/SVMCrossVal@96 83ab2cfd-5345-466c-8aeb-2b2739fb922d

Christoph Budziszewski authored on05/01/2009 18:25:16
Showing1 changed files
... ...
@@ -1,17 +1,17 @@
1 1
 % function [decodePerformance rawTimecourse ] = calculateDecodePerformance(des,timeLineStart, timeLineEnd, decodeDuration, svmargs, conditionList, sessionList, voxelList, classList, labelMap,normalize)
2
-function outputStruct = calculateDecodePerformance(inputStruct)
2
+function outputStruct = calculateDecodePerformance(inputStruct,SubjectID)
3 3
 
4 4
 addpath 'libsvm-mat-2.88-1';
5 5
 
6 6
 outputStruct = struct;
7 7
 
8
-des             = inputStruct.des;
8
+des             = inputStruct.(SubjectID).des;
9 9
 timeLineStart   = inputStruct.frameShiftStart;
10 10
 timeLineEnd     = inputStruct.frameShiftEnd;
11 11
 decodeDuration  = inputStruct.decodeDuration;
12 12
 svmargs         = inputStruct.svmargs;
13 13
 sessionList     = inputStruct.sessionList;
14
-voxelList       = inputStruct.voxelList;
14
+voxelList       = inputStruct.(SubjectID).voxelList;
15 15
 % classList       = inputStruct.classList;
16 16
 % labelMap        = inputStruct.labelMap;
17 17
 smoothed       = inputStruct.smoothed;
Browse code

killed unnecessary assignin raw axis fixed value added smooth option

git-svn-id: https://svn.discofish.de/MATLAB/spmtoolbox/SVMCrossVal@95 83ab2cfd-5345-466c-8aeb-2b2739fb922d

Christoph Budziszewski authored on18/12/2008 15:29:41
Showing1 changed files
... ...
@@ -14,7 +14,7 @@ sessionList     = inputStruct.sessionList;
14 14
 voxelList       = inputStruct.voxelList;
15 15
 % classList       = inputStruct.classList;
16 16
 % labelMap        = inputStruct.labelMap;
17
-% normalize       = inputStruct.normalize;
17
+smoothed       = inputStruct.smoothed;
18 18
 globalStart     = inputStruct.psthStart;
19 19
 globalEnd       = inputStruct.psthEnd;
20 20
 baselineStart   = inputStruct.baselineStart;
... ...
@@ -39,7 +39,7 @@ maxPerformance = -inf;
39 39
         
40 40
         
41 41
         for voxel = 1:size(voxelList,1)  % [[x;x],[y;y],[z;z]]
42
-                extr  = calculateImageData(voxelList(voxel,:),des);
42
+                extr  = calculateImageData(voxelList(voxel,:),des,smoothed);
43 43
                 rawdata=cell2mat({extr.mean}); % Raw Data
44 44
                 pst{voxel}  = calculatePST(des,globalStart,baselineStart,baselineEnd,globalEnd,eventList,rawdata,sessionList);
45 45
         end
... ...
@@ -83,7 +83,7 @@ maxPerformance = -inf;
83 83
 end
84 84
 
85 85
 
86
-function extr = calculateImageData(voxelList,des)
86
+function extr = calculateImageData(voxelList,des,smoothed)
87 87
 
88 88
 dtype='PSTH';
89 89
 
... ...
@@ -93,9 +93,15 @@ switch dtype
93 93
     case 'betas'
94 94
         V=des.Vbeta;
95 95
 end;
96
-%   for z=1:length(V) % Change Drive Letter!
97
-%       V(z).fname(1)='E';
98
-%   end;
96
+
97
+if (~smoothed)
98
+  for z=1:length(V) % Change Drive Letter!\
99
+      % D:....SUBJECTID\session\swfanders...
100
+      % D:....SUBJECTID\session\wfanders...
101
+      tmp = findstr(filesep,V(z).fname);
102
+      V(z).fname = strcat(V(z).fname(1:tmp(end)),V(z).fname(tmp(end)+2:end));
103
+  end;
104
+end
99 105
 
100 106
 % rad = 0; % one voxel
101 107
 % opt = 1; % xyz coordinates [mm]
Browse code

SVMCrossVal toolbox init

git-svn-id: https://svn.discofish.de/MATLAB/spmtoolbox/SVMCrossVal@90 83ab2cfd-5345-466c-8aeb-2b2739fb922d

Christoph Budziszewski authored on17/12/2008 13:45:29
Showing1 changed files
1 1
new file mode 100644
... ...
@@ -0,0 +1,278 @@
1
+% function [decodePerformance rawTimecourse ] = calculateDecodePerformance(des,timeLineStart, timeLineEnd, decodeDuration, svmargs, conditionList, sessionList, voxelList, classList, labelMap,normalize)
2
+function outputStruct = calculateDecodePerformance(inputStruct)
3
+
4
+addpath 'libsvm-mat-2.88-1';
5
+
6
+outputStruct = struct;
7
+
8
+des             = inputStruct.des;
9
+timeLineStart   = inputStruct.frameShiftStart;
10
+timeLineEnd     = inputStruct.frameShiftEnd;
11
+decodeDuration  = inputStruct.decodeDuration;
12
+svmargs         = inputStruct.svmargs;
13
+sessionList     = inputStruct.sessionList;
14
+voxelList       = inputStruct.voxelList;
15
+% classList       = inputStruct.classList;
16
+% labelMap        = inputStruct.labelMap;
17
+% normalize       = inputStruct.normalize;
18
+globalStart     = inputStruct.psthStart;
19
+globalEnd       = inputStruct.psthEnd;
20
+baselineStart   = inputStruct.baselineStart;
21
+baselineEnd     = inputStruct.baselineEnd;
22
+eventList       = inputStruct.eventList;
23
+
24
+
25
+minPerformance = inf;
26
+maxPerformance = -inf;
27
+
28
+
29
+        
30
+        %Pro Voxel PSTH TIMELINE berechnen.
31
+        %   timeshift mit pst-timeline durchf�hren.
32
+        % psth-timeline -25 bis +15 zu RES Onset.
33
+        
34
+%         eventList       = [9,11,13;10,12,14];
35
+%         globalStart     = -25;
36
+%         globalEnd       = 15;
37
+%         baselineStart   = -22;
38
+%         baselineEnd     = -20;
39
+        
40
+        
41
+        for voxel = 1:size(voxelList,1)  % [[x;x],[y;y],[z;z]]
42
+                extr  = calculateImageData(voxelList(voxel,:),des);
43
+                rawdata=cell2mat({extr.mean}); % Raw Data
44
+                pst{voxel}  = calculatePST(des,globalStart,baselineStart,baselineEnd,globalEnd,eventList,rawdata,sessionList);
45
+        end
46
+
47
+        decodePerformance = [];
48
+
49
+        for timeShift   = timeLineStart:1:timeLineEnd
50
+            frameStart  = floor(-globalStart+1+timeShift - 0.5*decodeDuration);
51
+            frameEnd    = min(ceil(frameStart+decodeDuration + 0.5*decodeDuration),-globalStart+globalEnd);
52
+            
53
+            tmp =[];
54
+            anyvoxel = 1;
55
+            for label = 1:size(pst{1,anyvoxel},2) 
56
+                for dp = 1:size(pst{1,anyvoxel}{1,label},1) % data point
57
+                row = label;
58
+                    for voxel = 1:size(pst,2)
59
+                        row = [row, pst{1,voxel}{1,label}(dp,frameStart:frameEnd)]; % label,value,value
60
+                    end
61
+                tmp  = [tmp; row];
62
+                end
63
+            end 
64
+        
65
+            svmdata      = tmp(:,2:size(tmp,2));
66
+            svmlabel     = tmp(:,1);
67
+            performance  = svmtrain(svmlabel, svmdata, svmargs);
68
+
69
+            minPerformance = min(minPerformance,performance);
70
+            maxPerformance = max(maxPerformance,performance);
71
+
72
+            decodePerformance = [decodePerformance; performance];
73
+        end
74
+        
75
+        outputStruct.decodePerformance  = decodePerformance;
76
+        outputStruct.svmdata            = svmdata;
77
+        outputStruct.svmlabel           = svmlabel;
78
+        outputStruct.rawTimeCourse      = pst;
79
+        outputStruct.minPerformance     = minPerformance;
80
+        outputStruct.maxPerformance     = maxPerformance;
81
+
82
+% display(sprintf('Min CrossVal Accuracy: %g%% \t Max CrossVal Accuracy: %g%%',minPerformance,maxPerformance));
83
+end
84
+
85
+
86
+function extr = calculateImageData(voxelList,des)
87
+
88
+dtype='PSTH';
89
+
90
+switch dtype 
91
+    case 'PSTH'
92
+        V=des.xY.VY;
93
+    case 'betas'
94
+        V=des.Vbeta;
95
+end;
96
+%   for z=1:length(V) % Change Drive Letter!
97
+%       V(z).fname(1)='E';
98
+%   end;
99
+
100
+% rad = 0; % one voxel
101
+% opt = 1; % xyz coordinates [mm]
102
+
103
+
104
+vox = voxelList;
105
+nRoi = size(vox,1);
106
+
107
+nImg = numel(V);
108
+
109
+for k=1:nImg
110
+	extr(k) = struct(...
111
+        'val',   repmat(NaN, [1 nRoi]),...
112
+		'mean',  repmat(NaN, [1 nRoi]),...
113
+		'sum',   repmat(NaN, [1 nRoi]),...
114
+		'nvx',   repmat(NaN, [1 nRoi]),...
115
+		'posmm', repmat(NaN, [3 nRoi]),...
116
+		'posvx', repmat(NaN, [3 nRoi]));
117
+
118
+    roicenter = round(inv(V(k).mat)*[vox, ones(nRoi,1)]');
119
+
120
+	for l = 1:nRoi
121
+
122
+%         if rad==0
123
+            x = roicenter(1,l);
124
+            y = roicenter(2,l);
125
+            z = roicenter(3,l);
126
+%         else
127
+%             tmp = spm_imatrix(V(k).mat);
128
+%             vdim = tmp(7:9);
129
+%             vxrad = ceil((rad*ones(1,3))./(ones(nRoi,1)*vdim))';
130
+%             [x y z] = ndgrid(-vxrad(1,l):sign(vdim(1)):vxrad(1,l), ...
131
+%                       -vxrad(2,l):sign(vdim(2)):vxrad(2,l), ...
132
+%                       -vxrad(3,l):sign(vdim(3)):vxrad(3,l));
133
+%             sel = (x./vxrad(1,l)).^2 + (y./vxrad(2,l)).^2 + ...
134
+%                   (z./vxrad(3,l)).^2 <= 1;
135
+%             x = roicenter(1,l)+x(sel(:));
136
+%             y = roicenter(2,l)+y(sel(:));
137
+%             z = roicenter(3,l)+z(sel(:));
138
+%         end;
139
+		dat                 = spm_sample_vol(V(k), x, y, z,0);
140
+		[maxv maxi]         = max(dat);
141
+		tmp                 = V(k).mat*[x(maxi); y(maxi); z(maxi);1]; % Max Pos
142
+		extr(k).val(l)      = maxv;
143
+		extr(k).sum(l)      = sum(dat);
144
+		extr(k).mean(l)     = nanmean(dat);
145
+        extr(k).nvx(l)      = numel(dat);
146
+		extr(k).posmm(:,l)  = tmp(1:3);
147
+		extr(k).posvx(:,l)  = [x(maxi); y(maxi); z(maxi)]; % Max Pos
148
+	end;
149
+
150
+end;
151
+end
152
+
153
+% disp(sprintf('Extracted at %.1f %.1f %.1f [xyz(mm)], average of %i voxel(s) [%.1fmm radius Sphere]',vox,length(x),rad));
154
+
155
+function pst = calculatePST(des,globalStart,baselineStart,baselineEnd,globalEnd,eventList,data,sessionList)
156
+    bstart          = baselineStart;
157
+    bend            = baselineEnd;
158
+    edur            = 12;
159
+    pre             =  globalStart;
160
+    post            =  globalEnd;
161
+    res             = 1;
162
+
163
+    normz           = 'file';
164
+    pm              = 0;
165
+
166
+    lsess           = getNumberOfScans(des);
167
+    nSessions       = getNumberOfSessions(des);
168
+    tr              = 2;
169
+
170
+    [evntrow evntcol]=size(eventList);
171
+    
172
+
173
+    hsec=str2num(des.xsDes.High_pass_Filter(8:end-3)); % Highpass filter [sec] from SPM.mat
174
+
175
+    if strcmp(des.xBF.UNITS,'secs')
176
+        unitsecs=1;
177
+    end;
178
+
179
+    nScansPerSession=getNumberOfScans(des);
180
+    %stime=[0:tr:max(nScansPerSession)*tr+post-tr]; % Stimulus time for raw data plot
181
+    stime=0:tr:max(nScansPerSession)*tr+round(post/tr)*tr-tr; % Stimulus time for raw data plot
182
+
183
+
184
+
185
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
186
+    % RUN
187
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
188
+
189
+
190
+    % Digital Highpass
191
+    Rp=0.5;
192
+    Rs=20;
193
+    NO=1;
194
+    Wp=1/((1/2/tr)/(1/hsec));
195
+    [B, A] = ellip(NO,Rp,Rs,Wp,'high');
196
+
197
+    sdata(1:max(nScansPerSession)+round(post/tr),1:nSessions)=nan; % Open Data Matrix
198
+    for z=1:nSessions % Fill Data Matrix sessionwise
199
+        sdata(1:nScansPerSession(z),z)=data(sum(nScansPerSession(1:z))-nScansPerSession(z)+1:sum(nScansPerSession(1:z)))';
200
+    end;
201
+%         usdata=sdata; % Keep unfiltered data
202
+
203
+    sdatamean=nanmean(nanmean(sdata(:,:)));
204
+    for z=1:nSessions
205
+%             X(:,z)=[1:1:max(nScansPerSession)]'; % #Volume
206
+        sdata(1:nScansPerSession(z),z)=filtfilt(B,A,sdata(1:nScansPerSession(z),z)); %Filter Data (Highpass)
207
+    end;
208
+    sdata=sdata+sdatamean;
209
+
210
+
211
+    %%%%Parametric Modulation Modus%%%%
212
+    if pm %Find Parameters for Event of Interest
213
+        [imods modss mods erow evntrow eventList] = getParametricMappingEvents(eventList,evntrow,des,pmf);
214
+    end;
215
+    %%%%PM%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
216
+
217
+
218
+    for zr=1:evntrow
219
+        n{zr}=0;
220
+        nn{zr}=0; 
221
+        nnn{zr}=0;
222
+        sstart{zr}=1;
223
+    end;
224
+
225
+
226
+    sesst0=0; 
227
+    for sessionID=sessionList
228
+        if sessionID>1
229
+            sesst0(sessionID)=sum(lsess(1:sessionID-1))*tr;  
230
+        end;
231
+        for zr=1:evntrow  %LABEL NUMBER, EVENT GROUP
232
+            sstart{zr}=n{zr}+1;
233
+            for ze=1:evntcol % EVENT INDEX in EventList
234
+                if ze==1 || (ze>1 && eventList(zr,ze)~=eventList(zr,ze-1))
235
+                    for zz=1:length(des.Sess(sessionID).U(eventList(zr,ze)).ons) % EVENT REPETITION NUMBER
236
+                        if ~unitsecs
237
+                            des.Sess(sessionID).U(eventList(zr,ze)).ons(zz)=(des.Sess(sessionID).U(eventList(zr,ze)).ons(zz)-1)*tr;
238
+                            des.Sess(sessionID).U(eventList(zr,ze)).dur(zz)=(des.Sess(sessionID).U(eventList(zr,ze)).dur(zz)-1)*tr;
239
+                        end;
240
+
241
+                        nnn{zr}=nnn{zr}+1; % INFO for rawdataplot start
242
+                        if des.Sess(sessionID).U(eventList(zr,ze)).dur(zz)<edur
243
+                            mev{zr}(nnn{zr},1:2)=[des.Sess(sessionID).U(eventList(zr,ze)).ons(zz)+sesst0(sessionID) edur]; % modeled event [onset length]
244
+                        else
245
+                            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)];
246
+                        end; % INFO for rawdataplot end
247
+
248
+                        n{zr}=n{zr}+1;
249
+                        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');
250
+                        if strcmp(normz,'epoc')
251
+                            bline=nanmean(pst{zr}(n{zr},round(-pre/res+(bstart)/res+1):round(-pre/res+(bend)/res+1)));
252
+                            if isnan(bline)
253
+                                pst{zr}(n{zr},1:-pre/res+post/res+1)=nan;
254
+                            else
255
+%                                     nn{zr}=nn{zr}+1;
256
+                                pst{zr}(n{zr},:)=(pst{zr}(n{zr},:)-bline)/bline*100; % 'epoch-based' normalization
257
+                            end;
258
+                        end;
259
+                    end;
260
+                end;
261
+            end;
262
+            if ~strcmp(normz,'epoc')
263
+                bline(zr)=nanmean(nanmean(pst{zr}(sstart{zr}:n{zr},-pre/res+(bstart)/res+1:-pre/res+(bend)/res+1)));
264
+                bstd(zr)=nanmean(nanstd(pst{zr}(sstart{zr}:n{zr},-pre/res+(bstart)/res+1:-pre/res+(bend)/res+1)));
265
+                nn{zr}=n{zr};
266
+            end;
267
+        end;
268
+        if strcmp(normz,'filz')
269
+            for zr=1:evntrow
270
+                pst{zr}(sstart{zr}:n{zr},:)=(pst{zr}(sstart{zr}:n{zr},:)-mean(bline))/mean(bstd); % session-based z-score normalization
271
+            end;
272
+        elseif strcmp(normz,'file')
273
+            for zr=1:evntrow
274
+                pst{zr}(sstart{zr}:n{zr},:)=(pst{zr}(sstart{zr}:n{zr},:)-mean(bline))/mean(bline)*100; % session-based normalization
275
+            end;
276
+        end;
277
+    end;
278
+end