refactored to arg-structs! (untested)
Christoph Budziszewski

Christoph Budziszewski commited on 2009-01-29 16:25:27
Zeige 9 geänderte Dateien mit 223 Einfügungen und 102 Löschungen.


git-svn-id: https://svn.discofish.de/MATLAB/spmtoolbox/SVMCrossVal@122 83ab2cfd-5345-466c-8aeb-2b2739fb922d
... ...
@@ -1,14 +1,14 @@
1
-function timePointMatrix = buildTimePointMatrix(argStruct)
1
+function timePointMatrix = buildTimePointMatrix(timeline,argStruct)
2 2
 
3
-pst             = argStruct.pst;
4 3
 
5
-timeLineStart   = argStruct.timeLineStart;
6
-timeLineEnd     = argStruct.timeLineEnd;
7
-globalStart     = argStruct.globalStart;
8
-globalEnd       = argStruct.globalEnd;
9
-decodeDuration  = argStruct.decodeDuration;
10
-eventList       = argStruct.eventList;
4
+timeLineStart   = timeline.frameShiftStart;
5
+timeLineEnd     = timeline.frameShiftEnd;
6
+globalStart     = timeline.psthStart;
7
+globalEnd       = timeline.psthEnd;
8
+decodeDuration  = timeline.decodeDuration;
11 9
 
10
+pst             = argStruct.pst;
11
+eventList       = argStruct.eventList;
12 12
 labelMap        = argStruct.labelMap;
13 13
 
14 14
 timePointMatrix = {};
... ...
@@ -1,5 +1,5 @@
1 1
 % function [decodePerformance rawTimecourse ] = calculateDecodePerformance(des,timeLineStart, timeLineEnd, decodeDuration, svmargs, conditionList, sessionList, voxelList, classList, labelMap,normalize)
2
-function outputStruct = calculateDecodePerformance(timeline,inputStruct,SubjectID)
2
+function outputStruct = calculateDecodePerformance(timeline,inputStruct,subjectParams)
3 3
 
4 4
 global CROSSVAL_METHOD_DEF;
5 5
 
... ...
@@ -10,13 +11,16 @@ METHOD              = inputStruct.CROSSVAL_METHOD;
10 11
 
11 12
 RANDOMIZE_DATAPOINTS = inputStruct.RANDOMIZE;
12 13
 
14
+% SubjectID       = subjectParams.SubjectID;
15
+% namehelper      = subjectParams.namehelper;
16
+voxelList       = subjectParams.voxelList;
17
+des             = subjectParams.des;
18
+
13 19
 outputStruct = struct;
14 20
 
15
-namehelper      = strcat('s',SubjectID);
16
-des             = inputStruct.(namehelper).des;
17 21
 svmargs         = inputStruct.svmargs;
18 22
 sessionList     = inputStruct.sessionList;
19
-voxelList       = inputStruct.(namehelper).voxelList;
23
+
20 24
 % classList       = inputStruct.classList;
21 25
 % labelMap        = inputStruct.labelMap;
22 26
 smoothed        = inputStruct.smoothed;
... ...
@@ -24,33 +28,45 @@ eventList       = inputStruct.eventList;
24 28
 
25 29
 timeLineStart   = timeline.frameShiftStart;
26 30
 timeLineEnd     = timeline.frameShiftEnd;
27
-decodeDuration  = timeline.decodeDuration;
28
-globalStart     = timeline.psthStart;
29
-globalEnd       = timeline.psthEnd;
30
-baselineStart   = timeline.baselineStart;
31
-baselineEnd     = timeline.baselineEnd;
31
+% decodeDuration  = timeline.decodeDuration;
32
+% globalStart     = timeline.psthStart;
33
+% globalEnd       = timeline.psthEnd;
34
+% baselineStart   = timeline.baselineStart;
35
+% baselineEnd     = timeline.baselineEnd;
32 36
 
33 37
 
34 38
 minPerformance = inf;
35 39
 maxPerformance = -inf;
36 40
 
37
-        %% ERSETZEN DURCH ROI-IMAGE!
38
-        for voxel = 1:size(voxelList,1)  % [[x;x],[y;y],[z;z]]
39
-                extr        = calculateImageData(voxelList(voxel,:),des,smoothed); 
40
-                rawdata     = cell2mat({extr.mean}); % Raw Data
41
-                pst{voxel}  = calculatePST(des,globalStart,baselineStart,baselineEnd,globalEnd,eventList,rawdata,sessionList);
41
+imageFiles = getImageFileList(des,~smoothed);
42
+
43
+
44
+extr = calculateImageData(imageFiles,voxelList);
45
+rawdata = extr.dat;
46
+
47
+nVoxel = size(voxelList,1);
48
+
49
+calculatePstOpts = struct;
50
+calculatePstOpts.des = des;
51
+calculatePstOpts.eventList = eventList;
52
+calculatePstOpts.sessionList = sessionList;
53
+
54
+for iVoxel = 1:nVoxel
55
+    pst{nVoxel} = calculatePST(timeline,calculatePstOpts,rawdata(iVoxel));
42 56
 end
43 57
 
58
+%         for voxel = 1:size(voxelList,1)  % [[x;x],[y;y],[z;z]]
59
+%             extr        = calculateImageData(imageFiles,voxelList(voxel,:));
60
+%             rawdata     = cell2mat({extr.mean}); % Raw Data
61
+%             pst{voxel}  = calculatePST(des,globalStart,baselineStart,baselineEnd,globalEnd,eventList,rawdata,sessionList);
62
+%         end
63
+
44 64
 timePointArgs.pst = pst;
45
-        timePointArgs.timeLineStart = timeLineStart;
46
-        timePointArgs.timeLineEnd   = timeLineEnd;
47
-        timePointArgs.globalStart   = globalStart;
48
-        timePointArgs.globalEnd     = globalEnd;
49
-        timePointArgs.decodeDuration= decodeDuration;
65
+
50 66
 timePointArgs.labelMap      = inputStruct.labelMap;
51 67
 timePointArgs.eventList     = eventList;
52 68
 
53
-        timePointMatrix = buildTimePointMatrix(timePointArgs);
69
+timePointMatrix = buildTimePointMatrix(timeline,timePointArgs);
54 70
 
55 71
 decodePerformance = [];
56 72
 for index = 1:timeLineEnd-timeLineStart+1
... ...
@@ -1,79 +1,102 @@
1
-function extr = calculateImageData(voxelList,des,smoothed)
1
+% filenameList as in SPM.xY.VY
2
+% voxelList in [x y z], 1 coordinate per row
2 3
 
3
-dtype='PSTH';
4
+function extr = calculateImageData(filenameList, voxelList)
5
+%function extr = calculateImageData(voxelList,des,smoothed)
6
+% global USE_DRIVE_CHECK_HACK;
4 7
 
5
-switch dtype 
6
-    case 'PSTH'
7
-        V=des.xY.VY;
8
-    case 'betas'
9
-        V=des.Vbeta;
10
-end;
8
+% dtype='PSTH';
11 9
 
12
-if V(1).fname(1)~='D'
13
-     for z=1:length(V) % Change Drive Letter!\
14
-      	V(z).fname(1) = 'D';
15
-     end; 
16
-end
10
+% switch dtype 
11
+%     case 'PSTH'
12
+%         V=des.xY.VY;
13
+%     case 'betas'
14
+%         V=des.Vbeta;
15
+% end;
17 16
 
18
-if (~smoothed)
19
-  for z=1:length(V) % Change Drive Letter!\
20
-      % D:....SUBJECTID\session\swfanders...
21
-      % D:....SUBJECTID\session\wfanders...
22
-      tmp = findstr(filesep,V(z).fname);
23
-      V(z).fname = strcat(V(z).fname(1:tmp(end)),V(z).fname(tmp(end)+2:end));
24
-  end;
25
-end
17
+% if USE_DRIVE_CHECK_HACK
18
+% if V(1).fname(1)~='D'
19
+%      for z=1:length(V) % Change Drive Letter - HACK!
20
+%       	V(z).fname(1) = 'D';
21
+%      end; 
22
+% end
23
+% end
26 24
 
27
-% rad = 0; % one voxel
28
-% opt = 1; % xyz coordinates [mm]
25
+% if (~smoothed)
26
+%   for z=1:length(V) % Change smoothed Filename - HACK!
27
+%       % D:....SUBJECTID\session\swfanders...
28
+%       % D:....SUBJECTID\session\wfanders...
29
+%       tmp = findstr(filesep,V(z).fname);
30
+%       V(z).fname = strcat(V(z).fname(1:tmp(end)),V(z).fname(tmp(end)+2:end));
31
+%   end;
32
+% end
29 33
 
34
+V = filenameList;
30 35
 
31 36
 vox = voxelList;
32
-nRoi = size(vox,1);
37
+nVoxel = size(vox,1);
38
+nImage = numel(V);
39
+
40
+for kImage=1:nImage
41
+	extr(kImage) = struct(...
42
+        'val',   repmat(NaN, [1 nVoxel]),...
43
+		'mean',  repmat(NaN, [1 nVoxel]),...
44
+		'sum',   repmat(NaN, [1 nVoxel]),...
45
+		'nvx',   repmat(NaN, [1 nVoxel]),...
46
+		'posmm', repmat(NaN, [3 nVoxel]),...
47
+		'posvx', repmat(NaN, [3 nVoxel]));
33 48
 
34
-nImg = numel(V);
49
+    switch TRANSFORM_METHOD
50
+        case 'single'
51
+            roicenter = round(inv(V(kImage).mat)*[vox, ones(nVoxel,1)]');
52
+            x = roicenter(1,:);
53
+            y = roicenter(2,:);
54
+            z = roicenter(3,:);
55
+        case 'image'
56
+             x = []; y = []; z = [];
57
+            for iRoiFile = 1:nRoiFiles 
58
+                [x1 y1] = ndgrid(1:V(k).dim(1),1:V(k).dim(2));
59
+                for p = 1:V(k).dim(3) % resample mask Vm(iRoiFile) in space of V(k)
60
+                    B = spm_matrix([0 0 -p 0 0 0 1 1 1]);
61
+                    M = inv(B*inv(V(k).mat)*Vm(iRoiFile).mat);
62
+                    msk = find(spm_slice_vol(Vm(iRoiFile),M,V(k).dim(1:2),0));
63
+                    if ~isempty(msk)
64
+                        z1 = p*ones(size(msk(:)));
65
+                        x = [x; x1(msk(:))];
66
+                        y = [y; y1(msk(:))];
67
+                        z = [z; z1];
68
+                    end;
69
+                end;
35 70
 
36
-for k=1:nImg
37
-	extr(k) = struct(...
38
-        'val',   repmat(NaN, [1 nRoi]),...
39
-		'mean',  repmat(NaN, [1 nRoi]),...
40
-		'sum',   repmat(NaN, [1 nRoi]),...
41
-		'nvx',   repmat(NaN, [1 nRoi]),...
42
-		'posmm', repmat(NaN, [3 nRoi]),...
43
-		'posvx', repmat(NaN, [3 nRoi]));
71
+            end
72
+    end
44 73
     
45
-    roicenter = round(inv(V(k).mat)*[vox, ones(nRoi,1)]');
46 74
 
47
-	for l = 1:nRoi
75
+	for iVoxel = 1:nVoxel
48 76
 
49 77
 %         if rad==0
50
-            x = roicenter(1,l);
51
-            y = roicenter(2,l);
52
-            z = roicenter(3,l);
78
+            x = roicenter(1,iVoxel);
79
+            y = roicenter(2,iVoxel);
80
+            z = roicenter(3,iVoxel);
53 81
 %         else
54
-%             tmp = spm_imatrix(V(k).mat);
82
+%             tmp = spm_imatrix(V(kImage).mat);
55 83
 %             vdim = tmp(7:9);
56
-%             vxrad = ceil((rad*ones(1,3))./(ones(nRoi,1)*vdim))';
57
-%             [x y z] = ndgrid(-vxrad(1,l):sign(vdim(1)):vxrad(1,l), ...
58
-%                       -vxrad(2,l):sign(vdim(2)):vxrad(2,l), ...
59
-%                       -vxrad(3,l):sign(vdim(3)):vxrad(3,l));
60
-%             sel = (x./vxrad(1,l)).^2 + (y./vxrad(2,l)).^2 + ...
61
-%                   (z./vxrad(3,l)).^2 <= 1;
62
-%             x = roicenter(1,l)+x(sel(:));
63
-%             y = roicenter(2,l)+y(sel(:));
64
-%             z = roicenter(3,l)+z(sel(:));
84
+%             vxrad = ceil((rad*ones(1,3))./(ones(nVoxel,1)*vdim))';
85
+%             [x y z] = ndgrid(-vxrad(1,iVoxel):sign(vdim(1)):vxrad(1,iVoxel), ...
86
+%                       -vxrad(2,iVoxel):sign(vdim(2)):vxrad(2,iVoxel), ...
87
+%                       -vxrad(3,iVoxel):sign(vdim(3)):vxrad(3,iVoxel));
88
+%             sel = (x./vxrad(1,iVoxel)).^2 + (y./vxrad(2,iVoxel)).^2 + ...
89
+%                   (z./vxrad(3,iVoxel)).^2 <= 1;
90
+%             x = roicenter(1,iVoxel)+x(sel(:));
91
+%             y = roicenter(2,iVoxel)+y(sel(:));
92
+%             z = roicenter(3,iVoxel)+z(sel(:));
65 93
 %         end;
66 94
 
67 95
 
68
-		dat                 = spm_sample_vol(V(k), x, y, z,0);
69
-		[maxv maxi]         = max(dat);
70
-		tmp                 = V(k).mat*[x(maxi); y(maxi); z(maxi);1]; % Max Pos
71
-		extr(k).val(l)      = maxv;
72
-		extr(k).sum(l)      = sum(dat);
73
-		extr(k).mean(l)     = nanmean(dat);
74
-        extr(k).nvx(l)      = numel(dat);
75
-		extr(k).posmm(:,l)  = tmp(1:3);
76
-		extr(k).posvx(:,l)  = [x(maxi); y(maxi); z(maxi)]; % Max Pos
96
+		dat = spm_sample_vol(V(kImage), x, y, z,0);
97
+        extr(kImage).dat(iVoxel)      = dat;
98
+		extr(kImage).mean(iVoxel)     = nanmean(dat);
99
+        extr(kImage).nvx(iVoxel)      = numel(dat);
77 100
 	end;
78 101
 
79 102
 end;
... ...
@@ -1,9 +1,14 @@
1
-function pst = calculatePST(des,globalStart,baselineStart,baselineEnd,globalEnd,eventList,data,sessionList)
2
-    bstart          = baselineStart;
3
-    bend            = baselineEnd;
1
+%function pst = calculatePST(des,globalStart,baselineStart,baselineEnd,globalEnd,eventList,data,sessionList)
2
+function pst = calculatePST(timeline,pstopts,data)
3
+    des         = pstopts.des;
4
+    eventList   = pstopts.eventList;
5
+    sessionList = pstopts.sessionList;
6
+
7
+    bstart          = timeline.baselineStart;
8
+    bend            = timeline.baselineEnd;
4 9
     edur            = 12;
5
-    pre             =  globalStart;
6
-    post            =  globalEnd;
10
+    pre             = timeline.psthStart;
11
+    post            = timeline.psthEnd;
7 12
     res             = 1;
8 13
 
9 14
     normz           = 'file';
... ...
@@ -1,6 +1,7 @@
1 1
 function classify(varargin)
2 2
 
3 3
 global CROSSVAL_METHOD_DEF;
4
+global VOXEL_SELECTION_MODE_DEF;
4 5
 
5 6
 switch nargin
6 7
     case 1
... ...
@@ -28,6 +29,7 @@ end
28 29
         
29 30
 %         calculateParams.CROSSVAL_METHOD_DEF = CROSSVAL_METHOD_DEF;
30 31
         calculateParams.CROSSVAL_METHOD      = CROSSVAL_METHOD_DEF.svmcrossval;
32
+        calculateParams.VOXEL_SELECTION_MODE = VOXEL_SELECTION_MODE_DEF.roiImage;
31 33
         calculateParams.PROJECT_BASE_PATH    = PROJECT_BASE_PATH;
32 34
         calculateParams.PROJECT_RESULT_PATH  = PROJECT_RESULT_PATH;
33 35
         
... ...
@@ -76,11 +78,13 @@ end
76 78
 
77 79
 %% subject loop
78 80
 function decode = claculateMultiSubjectDecodePerformance(timelineParams,calculateParams)
81
+global VOXEL_SELECTION_MODE_DEF;
82
+
79 83
 decode = struct;
80 84
 decode.decodePerformance = [];
81 85
 decode.rawTimeCourse     = [];
82 86
 
83
-for subjectCell = subjectSelection
87
+for subjectCell = calculateParams.subjectSelection
84 88
     SubjectID = cell2mat(subjectCell);
85 89
     namehelper = strcat('s',SubjectID); %Vars can not start with numbers.
86 90
 
... ...
@@ -88,16 +92,29 @@ for subjectCell = subjectSelection
88 92
     spm = load(fullfile(calculateParams.PROJECT_BASE_PATH,SubjectID,calculateParams.PROJECT_RESULT_PATH));
89 93
     display('... done.');
90 94
 
91
-    % calculate
92
-    calculateParams.(namehelper).des             = spm.SPM;
93
-    calculateParams.(namehelper).voxelList       = parseVoxelList(paramModel,SubjectID);
94
-    assignin('base','calculateParams',calculateParams);
95
+    % per subject params..
96
+    subjectParams = struct;
97
+    
98
+    subjectParams.des             = spm.SPM;
99
+
100
+    switch calculateParams.VOXEL_SELECTION_MODE 
101
+        case VOXEL_SELECTION_MODE_DEF.manualGui
102
+            subjectParams.voxelList       = parseVoxelList(paramModel,SubjectID);
103
+            
104
+        case VOXEL_SELECTION_MODE_DEF.roiImage
105
+            subjectParams.voxelList       = readRoiImage();
106
+    end
107
+    
108
+    
109
+    subjectParams.SubjectID       = SubjectID;
110
+    subjectParams.namehelper      = namehelper;
111
+%     assignin('base','calculateParams',calculateParams);
95 112
 
96 113
     display(sprintf('calculating cross-validation performance time-shift for Subject %s. Please Wait. ...',SubjectID));
97 114
     display('switching off all warnings');
98 115
     warning_state               = warning('off','all');
99 116
     display('calculating ...');
100
-    decode.(namehelper)         = calculateDecodePerformance(timelineParams,calculateParams,SubjectID);
117
+    decode.(namehelper)         = calculateDecodePerformance(timelineParams,calculateParams,subjectParams);
101 118
 
102 119
     display('... done');
103 120
     display('restoring warnings');
... ...
@@ -0,0 +1,24 @@
1
+function fileList = getImageFileList(des,use_unsmoothed_hack)
2
+
3
+global USE_DRIVE_CHECK_HACK;
4
+
5
+fileList  = des.xY.VY;
6
+nFiles = length(fileList);
7
+
8
+if USE_DRIVE_CHECK_HACK
9
+    if fileList(1).fname(1)~='D'
10
+        for z=1:nFiles % Change Drive Letter - HACK!
11
+            fileList(z).fname(1) = 'D';
12
+        end;
13
+    end
14
+end
15
+if use_unsmoothed_hack
16
+    for z=1:nFiles % Change smoothed Filename - HACK!
17
+        % D:....SUBJECTID\session\swfanders...
18
+        % D:....SUBJECTID\session\wfanders...
19
+        tmp = findstr(filesep,fileList(z).fname);
20
+        fileList(z).fname = strcat(fileList(z).fname(1:tmp(end)),fileList(z).fname(tmp(end)+2:end));
21
+    end;
22
+end
23
+
24
+end
0 25
\ No newline at end of file
... ...
@@ -1,5 +1,7 @@
1 1
 function main_UI(args)
2 2
 
3
+global CROSSVAL_METHOD_DEF;
4
+
3 5
 %     project = args.project;
4 6
 %     DEFAULT = project.default;
5 7
 
... ...
@@ -61,12 +63,17 @@ DEFAULT.svmoptstring    = '-s 0 -t 0 -v 6 -c 1';
61 63
                     'Position',[firstColumn firstRow 0.66*frameWidth controlElementHeight*6]);
62 64
     set(model.subjectSelector,'BackgroundColor','w');
63 65
     
64
-    createLabel(pSubject,[0.68*frameWidth firstRow*2  0.25*frameWidth controlElementHeight],'Smooth Data');
65 66
     model.chkSmoothed     = uicontrol(pSubject,'Style','checkbox','Position',[0.68*frameWidth firstRow  0.25*frameWidth controlElementHeight],'Value',DEFAULT.smoothed);
67
+    createLabel(pSubject,[0.75*frameWidth firstRow  0.25*frameWidth controlElementHeight],'Smooth Data?');
68
+
69
+    createLabel(pSubject,[0.68*frameWidth firstRow*4  0.25*frameWidth controlElementHeight],'Analysis Method');
66 70
     
67
-    createLabel(pSubject,[0.68*frameWidth firstRow*4  0.25*frameWidth controlElementHeight],'Crossvalidation');
68
-    model.txtMultisubject = createTextField(pSubject,[0.68*frameWidth firstRow*3  0.25*frameWidth controlElementHeight],DEFAULT.multisubject);
69
-    set(model.txtMultisubject,'enable','off');
71
+    model.analysisMethodSelector = uicontrol(pSubject,'Style','popupmenu',...
72
+        'Position',[0.68*frameWidth firstRow*3  0.25*frameWidth controlElementHeight],...
73
+        'String',{CROSSVAL_METHOD_DEF.svmcrossval, CROSSVAL_METHOD_DEF.classPerformance},...
74
+        'Value',1);
75
+    set(model.analysisMethodSelector,'BackgroundColor','w');
76
+%     set(model.analysisMethodSelector,'enable','off');
70 77
 
71 78
     % PSTH
72 79
     firstColumn  = 5.00;
... ...
@@ -0,0 +1,24 @@
1
+function voxelList = readRoiImage()
2
+
3
+Vm = spm_vol(spm_select([1 Inf],'image','Select ROI image'));
4
+nRoiFiles = size(Vm,2);
5
+
6
+% for iRoiFile = 1:nRoiFiles
7
+%     x = []; y = []; z = [];
8
+%     [x1 y1] = ndgrid(1:V(k).dim(1),1:V(k).dim(2));
9
+%     for p = 1:V(k).dim(3) % resample mask Vm(iRoiFile) in space of V(k)
10
+%         B = spm_matrix([0 0 -p 0 0 0 1 1 1]);
11
+%         M = inv(B*inv(V(k).mat)*Vm(iRoiFile).mat);
12
+%         msk = find(spm_slice_vol(Vm(iRoiFile),M,V(k).dim(1:2),0));
13
+%         if ~isempty(msk)
14
+%             z1 = p*ones(size(msk(:)));
15
+%             x = [x; x1(msk(:))];
16
+%             y = [y; y1(msk(:))];
17
+%             z = [z; z1];
18
+%         end;
19
+%     end;
20
+% 
21
+% end
22
+voxelList = Vm;
23
+
24
+end
0 25
\ No newline at end of file
... ...
@@ -1,13 +1,19 @@
1 1
 function spm_SVMCrossVal(varargin)
2 2
 
3 3
 %define global constants
4
-global CROSSVAL_METHOD_DEF;
4
+global USE_DRIVE_CHECK_HACK;
5
+USE_DRIVE_CHECK_HACK = 1;  %enables subroutine to check if image path starts with 'D'
5 6
 
7
+global CROSSVAL_METHOD_DEF;
6 8
 CROSSVAL_METHOD_DEF.svmcrossval       = 'svm crossval';
7 9
 CROSSVAL_METHOD_DEF.classPerformance  = 'svm class performance';
8
-CROSSVAL_METHOD_DEF.crossSubject      = 'svm cross subject testing';
10
+CROSSVAL_METHOD_DEF.crossSubject      = 'svm across subject testing';
9 11
 CROSSVAL_METHOD_DEF.somTraining       = 'som Training';
10 12
 
13
+global VOXEL_SELECTION_MODE_DEF;
14
+VOXEL_SELECTION_MODE_DEF.manualGui    = 'manually defined in GUI';
15
+VOXEL_SELECTION_MODE_DEF.roiImage     = 'use ROI image and popup image selector';
16
+
11 17
 switch nargin
12 18
 case 0
13 19
 % 	project_UI;
14 20