SVMCrossVal.git

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

2) addpath(fullfile(getTbxPath,'NIFTI_20090325'));

3) 

4) SVMCROSSVAL_DEBUG = 0;

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

6) 

7) subjectResultDir = subjectdata{1}.dir;
8) savePath = fullfile(subjectResultDir,'results','fbs',datestr(now,30));

9) mkdir(savePath);
10) 

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

20) radiusList = fbsargs.radius;

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

33) 
34) 

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

38) 
39) 

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

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

44) 

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

49)     tic;

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

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

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

54)     tic 

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

56)     elapsed{s}.calcExtr = toc;

57)     

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

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

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

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

72)         indexToCoordMap.put(iVoxel,coord);

73)         coordToIndexMap.put(a,iVoxel);

74)     end
75)     

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

79)     elapsed{s}.mapping = toc;

80)     

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

87)     tic

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

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

102)     end

103)     if isDebug

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

114)     elapsed{s}.psth = toc;

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

117)     clear extr;

118)     %run searchlight

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

120)     tic

121)     

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

123)     % for each timeslice

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

129)     

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

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

134)     

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

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

138)     

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

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

147) 

148) if isempty( fbsargs.timeline )

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

151) else 

152)     %preferred!

153) 	fbsTimeLine = fbsargs.timeline;

154) end
155) 

156) for timeShiftIdx = fbsTimeLine

157)         

158)     % center timepoint && relative shift

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

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

161) 

162)         for rIdx = 1:length(radiusList)

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

164)         end
165)         

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

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

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

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

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

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

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

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

188)                         for voxel = sphere

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

192)                     end
193)                 end

194)                 
195)                 if RANDOMIZE_DATAPOINTS
196)                     %                 [svmlabel svmdata]
197)                     rndindex  = randperm(length(svmlabel));
198)                     svmdata   = svmdata(rndindex,:);
199)                     svmlabel  = svmlabel(rndindex);
200)                 end
201)                 decode = svm_single_crossval(svmlabel,svmdata,svmopts);
202)                 % save the decode value to the corresponding coordinate
203)                 
204)                 coord = mapping.indexToCoordMap.get(iVoxel);
205)                 x = coord(1);
206)                 y = coord(2);
207)                 z = coord(3);
208)                 

209)                 if NORM_DECODE
210)                     img3D{rIdx}(x,y,z) = ((decode/100)-0.5)*2; % range [-1:1]
211)                 else
212)                     img3D{rIdx}(x,y,z) = (decode/100); % range [0:1]
213)                 end
214) 

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

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

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

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

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

224)         end

225)     end %for each timeslice

226)     display('rastering done');

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

229)     elapsed{s}.decode = toc;

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

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

235) end
236) 
237) end