SVMCrossVal.git

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

2) addpath('NIFTI_20090325');
3) 

4) DEBUG = 0;
5) NORM_DECODE = 0; % 1: normalize to [-1:1], 0: normalize to [0:1]

6) 

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

8) mkdir(savePath);
9) 

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

19) radiusList = fbsargs.radius;

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

32) 
33) 

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

37) 
38) 

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

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

43) 

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

48)     tic;

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

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

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

53)     tic 

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

55)     elapsed{s}.calcExtr = toc;

56)     

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

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

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

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

71)         indexToCoordMap.put(iVoxel,coord);

72)         coordToIndexMap.put(a,iVoxel);

73)     end
74)     

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

78)     elapsed{s}.mapping = toc;

79)     

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

86)     tic

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

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

101)     end

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

113)     elapsed{s}.psth = toc;

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

116)     clear extr;

117)     %run searchlight

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

119)     tic

120)     

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

122)     % for each timeslice

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

128)     

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

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

133)     

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

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

137)     

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

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

146) 

147) if isempty( fbsargs.timeline )

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

150) else 

151)     %preferred!

152) 	fbsTimeLine = fbsargs.timeline;

153) end
154) 

155) for timeShiftIdx = fbsTimeLine

156)         

157)     % center timepoint && relative shift

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

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

160) 

161)         for rIdx = 1:length(radiusList)

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

163)         end
164)         

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

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

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

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

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

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

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

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

187)                         for voxel = sphere

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

191)                     end
192)                 end

193)                 
194)                 if RANDOMIZE_DATAPOINTS
195)                     %                 [svmlabel svmdata]
196)                     rndindex  = randperm(length(svmlabel));
197)                     svmdata   = svmdata(rndindex,:);
198)                     svmlabel  = svmlabel(rndindex);
199)                 end
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)                 if NORM_DECODE
209)                     img3D{rIdx}(x,y,z) = ((decode/100)-0.5)*2; % range [-1:1]
210)                 else
211)                     img3D{rIdx}(x,y,z) = (decode/100); % range [0:1]
212)                 end
213) 

214)                 
215)             end %for each radius
216)         end %for each voxel
217) 

218)         for rIdx = 1:length(radiusList)
219)             radius = radiusList(rIdx);

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

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

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

223)         end

224)     end %for each timeslice

225)     display('rastering done');

226)     display(sprintf('result images saved to %s',savePath));
227)     

228)     elapsed{s}.decode = toc;

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

230) 
231)     %save memory!
232)     clear pst; 
233)     clear mask_image;

234) end
235) 
236) end