SVMCrossVal.git

1) function output =  runFBSImageMaskMode(header,subjectdata,fbsargs)

2) addpath('NIFTI_20090325');
3) 

4) savePath = fullfile('output',datestr(now,30));

5) mkdir(savePath);
6) 

7) args = header.args;
8) 
9) subjects = subjectdata;
10) nSubjects = numel(subjects);
11) sessionlist = args.sessionList;
12) pstopts = args.psthOpts;
13) pstopts.eventList = args.eventList;
14) pstopts.sessionList = sessionlist;
15) 

16) radiusList = fbsargs.radius;

17) 
18) timeline = header.timeline;
19) timeline.frameShiftStart = header.frameShift.frameShiftStart;
20) timeline.frameShiftEnd   = header.frameShift.frameShiftEnd;
21) timeline.decodeDuration  = header.frameShift.decodeDuration;
22) 
23) timePointArgs.labelMap      = LabelMap(header.classDef.labelCells,header.classDef.conditionCells);
24) timePointArgs.eventList     = header.classDef.eventMatrix;
25) 
26) timeLineStart  = timeline.frameShiftStart;
27) timeLineEnd    = timeline.frameShiftEnd;
28) 

29) 
30) 

31) RANDOMIZE_DATAPOINTS = header.svmrnd;
32) svmopts = fbsargs.svmopts;
33) 

34) 
35) 

36) disp(sprintf('batch processing  %g subjects.',nSubjects));
37) 
38) for s = 1:nSubjects

39)     elapsed{s} = struct; %measure the timing

40) 

41)     disp(sprintf('processing subject %s.',subjects{s}.name));
42)     % load image data
43)    
44)     disp('fetching volume definitions, please wait');

45)     tic;

46)     volumes = spm_vol(getImageFileList(subjects{s}.dir,sessionlist,args.mask));

47)     elapsed{s}.loadList = toc;
48)     

49)     disp('computing volume values, please wait');

50)     tic 

51)     [extr x y z] = calculateRoiImageData(volumes,subjects{s}.roiFile);

52)     elapsed{s}.calcExtr = toc;

53)     

54)     clear volumes; %save memory ??
55)     disp('cleared volumes');
56)     tic

57)     
58)     nVoxel = size(extr(1).dat,1);
59) 
60)     indexToCoordMap = java.util.HashMap;
61)     coordToIndexMap = java.util.HashMap;
62)     for iVoxel = 1:nVoxel

63)         coord = [x(iVoxel) y(iVoxel) z(iVoxel)];

64)         a = java.util.Vector(3,0);
65)         a.add(0,coord(1));
66)         a.add(1,coord(2));
67)         a.add(2,coord(3));

68)         indexToCoordMap.put(iVoxel,coord);

69)         coordToIndexMap.put(a,iVoxel);

70)     end
71)     

72)     mapping.indexToCoordMap = indexToCoordMap;
73)     mapping.coordToIndexMap = coordToIndexMap;
74) 

75)     elapsed{s}.mapping = toc;

76)     

77)     % calculate psth
78)     
79)     warn = warning('off','all');
80)     pstopts.des = subjects{s}.des;
81) 
82)     disp(sprintf('computing psth for %g voxel.',nVoxel));

83)     tic

84) 
85)     for iVoxel = 1:nVoxel
86)         rawdata = [];
87)         for iImage = 1:length(extr);
88)             tmp = extr(iImage);
89)             rawdata = [rawdata tmp.dat(iVoxel)];
90)         end
91) %         coordNameHelper = strcat('v',num2str(iVoxel));
92)         pst{iVoxel} = calculatePST(args.timeline,pstopts,rawdata); % do not store in subjectStruct, so we can clear it later.
93) 

94) %          p = iVoxel/nVoxel*100;
95) %          if(mod(iVoxel,floor(nVoxel/100))==0)
96) %              sprintf(' %g%%\t complete',p);
97) %          end

98)     end

99)     elapsed{s}.psth = toc;

100)     disp('psth done');
101)     warning(warn);

102)     clear extr;

103)     %run searchlight

104)     pause(0.001) % flush system event queue (respond to ctrl-c)

105)     tic

106)     

107)     display(sprintf('rastering %g coordinates',nVoxel));

108)     % for each timeslice

109)     globalStart     = timeline.psthStart;
110)     globalEnd       = timeline.psthEnd;
111)     decodeDuration  = timeline.decodeDuration;
112)     labelMap        = timePointArgs.labelMap;
113)     eventList       = pstopts.eventList;

114)     

115)     res             = timeline.tr*timeline.trFactor;
116)     

117)     % tmp = spm_imatrix(V(kImage).mat); %
118)     % vdim = tmp(7:9); % Voxel-Size

119)     

120)     mask_image = load_nii(subjects{s}.roiFile.fname);
121)     

122)     vdim = mask_image.hdr.dime.pixdim(2:4);
123)     if mask_image.hdr.dime.pixdim(1) == 1
124)         vdim = vdim .* [-1 1 1];
125)     end
126) 
127)     display(sprintf('starting timesliceing'));
128)     

129)     nSamplePoints = ((timeLineEnd-timeLineStart)/res) +1;

130) 

131) if isempty( fbsargs.timeline )

132) 	fbsTimeLine = 1:nSamplePoints;
133)     fbsTimeLine = fbsTimeLine +globalStart;

134) else 

135)     %preferred!

136) 	fbsTimeLine = fbsargs.timeline;

137) end
138) 

139) for timeShiftIdx = fbsTimeLine

140)         

141)     % center timepoint && relative shift

142)     frameStartIdx  = floor(-globalStart+1+timeShiftIdx - 0.5*decodeDuration);

143)     frameEndIdx    = min(ceil(frameStartIdx+decodeDuration + 0.5*decodeDuration),-globalStart+globalEnd);

144) 

145)         for r = radiusList
146)             img3D{r} = zeros(size(mask_image)); %output image prepare
147)         end
148)         

149)         for iVoxel = 1:nVoxel % linear structure avoids 3D-Loop.

150)             if (mod(iVoxel,100)== 0)

151)                 display(sprintf('Status: %03u / %03u Timepoints, %05u / %05u Coordinates',find(fbsTimeLine == timeShiftIdx),length(fbsTimeLine),iVoxel,nVoxel));

152)                 pause(0.001) %flush system event queue
153)             end

154)             for radius = radiusList
155)                 % get surrounding coordinate-IDs within radius
156)                 sphere = fbs_buildSphere(mapping,iVoxel,radius,vdim);
157)                 
158)                 %build svm inputmatrix
159)                 svmdata = [];
160)                 svmlabel = [];
161)                 anyvoxel = 1;
162)                 for pstConditionGroup = 1:size(pst{1,anyvoxel},2)
163)                     for dp = 1:size(pst{1,anyvoxel}{1,pstConditionGroup},1) % data point
164)                         svmlabel = [svmlabel; lm_getSVMLabel(labelMap,eventList(pstConditionGroup,1))];
165)                         row = [];
166)                         for voxel = 1:size(sphere,2)
167)                             row = [row, pst{1,voxel}{1,pstConditionGroup}(dp,frameStartIdx:frameEndIdx)]; % label,[value,value,...],[value,value,...]...
168)                         end
169)                         svmdata  = [svmdata; row];

170)                     end
171)                 end

172)                 
173)                 if RANDOMIZE_DATAPOINTS
174)                     %                 [svmlabel svmdata]
175)                     rndindex  = randperm(length(svmlabel));
176)                     svmdata   = svmdata(rndindex,:);
177)                     svmlabel  = svmlabel(rndindex);
178)                 end
179)                 
180)                 decode = svm_single_crossval(svmlabel,svmdata,svmopts);
181)                 % save the decode value to the corresponding coordinate
182)                 
183)                 coord = mapping.indexToCoordMap.get(iVoxel);
184)                 x = coord(1);
185)                 y = coord(2);
186)                 z = coord(3);
187)                 
188)                 img3D{radius}(x,y,z) = ((decode/100)-0.5)*2; % range [-1:1]
189)                 
190)             end %for each radius
191)         end %for each voxel
192) 
193)         for radius = radiusList
194)             nii = make_nii(img3D{radius},vdim,mask_image.hdr.hist.originator(1:3),16,...
195)                 sprintf('decode performance, time relative to onset: %g to %g sec',frameStartIdx,frameEndIdx));
196)             save_nii(nii,fullfile(savePath,sprintf('%s-r%g-t%+03g',subjects{s}.name,radius,timeShiftIdx)));
197)         end

198)     end %for each timeslice

199)     display('rastering done');

200)     display(sprintf('result images saved to %s',savePath));
201)     

202)     elapsed{s}.decode = toc;

203)     assignin('base','timing',elapsed);

204) 
205)     %save memory!
206)     clear pst; 
207)     clear mask_image;

208) end
209) 
210) end