SVMCrossVal.git

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

2) addpath('NIFTI_20090325');
3) 

4) DEBUG = 1;
5) 

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

7) mkdir(savePath);
8) 

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

18) radiusList = fbsargs.radius;

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

31) 
32) 

33) RANDOMIZE_DATAPOINTS = header.svmrnd;
34) svmopts = fbsargs.svmopts;
35) 

36) 
37) 

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

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

42) 

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

47)     tic;

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

49)     elapsed{s}.loadList = toc;
50)     

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

52)     tic 

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

54)     elapsed{s}.calcExtr = toc;

55)     

56)     clear volumes; %save memory ??
57)     disp('cleared volumes');
58)     tic

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

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

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

70)         indexToCoordMap.put(iVoxel,coord);

71)         coordToIndexMap.put(a,iVoxel);

72)     end
73)     

74)     mapping.indexToCoordMap = indexToCoordMap;
75)     mapping.coordToIndexMap = coordToIndexMap;
76) 

77)     elapsed{s}.mapping = toc;

78)     

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

85)     tic

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

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

100)     end

101)     if DEBUG == 1
102)         figure;
103)         hold on;
104)         for i = 1:size(pst,2)
105)             plot(mean(pst{:,i}{1}),'r-');
106)             plot(mean(pst{:,i}{2}),'b-');
107)         end
108)         hold off;
109)     end
110)     
111)     

112)     elapsed{s}.psth = toc;

113)     disp('psth done');
114)     warning(warn);

115)     clear extr;

116)     %run searchlight

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

118)     tic

119)     

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

121)     % for each timeslice

122)     globalStart     = timeline.psthStart;
123)     globalEnd       = timeline.psthEnd;
124)     decodeDuration  = timeline.decodeDuration;
125)     labelMap        = timePointArgs.labelMap;
126)     eventList       = pstopts.eventList;

127)     

128)     res             = timeline.tr*timeline.trFactor;
129)     

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

132)     

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

134)     tmp = strfind(mask_image.fileprefix,filesep);
135)     maskname = mask_image.fileprefix(tmp(end)+1:end);

136)     

137)     vdim = mask_image.hdr.dime.pixdim(2:4);
138)     if mask_image.hdr.dime.pixdim(1) == 1
139)         vdim = vdim .* [-1 1 1];
140)     end
141) 
142)     display(sprintf('starting timesliceing'));
143)     

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

145) 

146) if isempty( fbsargs.timeline )

147) 	fbsTimeLine = 1:nSamplePoints;
148)     fbsTimeLine = fbsTimeLine +globalStart;

149) else 

150)     %preferred!

151) 	fbsTimeLine = fbsargs.timeline;

152) end
153) 

154) for timeShiftIdx = fbsTimeLine

155)         

156)     % center timepoint && relative shift

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

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

159) 

160)         for rIdx = 1:length(radiusList)

161)             img3D{rIdx}(:,:,:) = zeros(size(mask_image)); %output image prepare

162)         end
163)         

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

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

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

167)                 pause(0.001) %flush system event queue
168)             end

169)             for rIdx  = 1: length(radiusList)
170)                 radius = radiusList(rIdx);

171)                 % get surrounding coordinate-IDs within radius
172)                 sphere = fbs_buildSphere(mapping,iVoxel,radius,vdim);

173) %                 for i = 1: length(sphere)
174) %                      sphere(i )
175) %                      mapping.indexToCoordMap.get(sphere(i))
176) %                 end

177)                 
178)                 %build svm inputmatrix
179)                 svmdata = [];
180)                 svmlabel = [];
181)                 anyvoxel = 1;
182)                 for pstConditionGroup = 1:size(pst{1,anyvoxel},2)
183)                     for dp = 1:size(pst{1,anyvoxel}{1,pstConditionGroup},1) % data point
184)                         svmlabel = [svmlabel; lm_getSVMLabel(labelMap,eventList(pstConditionGroup,1))];
185)                         row = [];

186)                         for voxel = sphere

187)                             row = [row, pst{1,voxel}{1,pstConditionGroup}(dp,frameStartIdx:frameEndIdx)]; % label,[value,value,...],[value,value,...]...
188)                         end
189)                         svmdata  = [svmdata; row];

190)                     end
191)                 end

192)                 
193)                 if RANDOMIZE_DATAPOINTS
194)                     %                 [svmlabel svmdata]
195)                     rndindex  = randperm(length(svmlabel));
196)                     svmdata   = svmdata(rndindex,:);
197)                     svmlabel  = svmlabel(rndindex);
198)                 end

199)                 decode = 0;

200)                 decode = svm_single_crossval(svmlabel,svmdata,svmopts);
201)                 % save the decode value to the corresponding coordinate
202)                 
203)                 coord = mapping.indexToCoordMap.get(iVoxel);
204)                 x = coord(1);
205)                 y = coord(2);
206)                 z = coord(3);
207)                 

208)                 img3D{rIdx}(x,y,z) = ((decode/100)-0.5)*2; % range [-1:1]

209)                 
210)             end %for each radius
211)         end %for each voxel
212) 

213)         for rIdx = 1:length(radiusList)
214)             radius = radiusList(rIdx);

215)             nii = make_nii(img3D{rIdx},vdim,mask_image.hdr.hist.originator(1:3),16,...

216)                 sprintf('decode performance, time relative to onset: %g to %g sec',frameStartIdx,frameEndIdx));

217)             save_nii(nii,fullfile(savePath,sprintf('%s-%s-r%g-t%+03g',subjects{s}.name,maskname,radius,timeShiftIdx)));

218)         end

219)     end %for each timeslice

220)     display('rastering done');

221)     display(sprintf('result images saved to %s',savePath));
222)     

223)     elapsed{s}.decode = toc;

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

225) 
226)     %save memory!
227)     clear pst; 
228)     clear mask_image;

229) end
230) 
231) end