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) radius = 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) RANDOMIZE_DATAPOINTS = header.svmrnd;
30) svmopts = fbsargs.svmopts;
31) 

32) 
33) 

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

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

38) 

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

43)     tic;

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

45)     elapsed{s}.loadList = toc;
46)     

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

48)     tic 

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

50)     elapsed{s}.calcExtr = toc;

51)     

52)     clear volumes; %save memory ??
53)     disp('cleared volumes');
54)     tic

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

61)         coord = [x(iVoxel) y(iVoxel) z(iVoxel)];
62)         a = java.util.Vector(3,0);
63)         a.add(0,coord(1));
64)         a.add(1,coord(2));
65)         a.add(2,coord(3));

66)         indexToCoordMap.put(iVoxel,coord);

67)         coordToIndexMap.put(a,iVoxel);

68)     end
69)     

70)     mapping.indexToCoordMap = indexToCoordMap;
71)     mapping.coordToIndexMap = coordToIndexMap;
72) 

73)     elapsed{s}.mapping = toc;

74)     

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

81)     tic

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

93)          if(mod(iVoxel,floor(nVoxel/100))==0)

94)              sprintf(' %g%%\t complete',p);
95)          end
96)     end

97)     elapsed{s}.psth = toc;

98)     disp('psth done');
99)     warning(warn);

100)     clear extr;

101)     %run searchlight

102)     tic

103)     
104)     display(sprintf('rastering %g coordinates with approx. %g mm radius',nVoxel,radius));
105)     % for each timeslice

106)     globalStart     = timeline.psthStart;
107)     globalEnd       = timeline.psthEnd;
108)     decodeDuration  = timeline.decodeDuration;
109)     labelMap        = timePointArgs.labelMap;
110)     eventList       = pstopts.eventList;

111)     
112)     % tmp = spm_imatrix(V(kImage).mat); %
113)     % vdim = tmp(7:9); % Voxel-Size

114)     
115) 

116)     

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

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

126)     for timeShift = 1:timeLineEnd-timeLineStart+1
127)     % center timepoint && relative shift
128)     frameStart  = floor(-globalStart+1+timeShift - 0.5*decodeDuration);
129)     frameEnd    = min(ceil(frameStart+decodeDuration + 0.5*decodeDuration),-globalStart+globalEnd);

130) 

131)         img3D = zeros(size(mask_image)); %output image prepare

132) 

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

134)             % get surrounding coordinate-IDs within radius
135)             sphere = fbs_buildSphere(mapping,iVoxel,radius,vdim);
136) 
137)             %build svm inputmatrix
138)             svmdata = [];
139)             svmlabel = [];
140)             anyvoxel = 1;
141)             for pstConditionGroup = 1:size(pst{1,anyvoxel},2)
142)                 for dp = 1:size(pst{1,anyvoxel}{1,pstConditionGroup},1) % data point
143)                     svmlabel = [svmlabel; lm_getSVMLabel(labelMap,eventList(pstConditionGroup,1))];
144)                     row = [];
145)                     for voxel = 1:size(sphere,2)
146)                         row = [row, pst{1,voxel}{1,pstConditionGroup}(dp,frameStart:frameEnd)]; % label,[value,value,...],[value,value,...]...
147)                     end
148)                     svmdata  = [svmdata; row];
149)                 end
150)             end

151) 
152)             if RANDOMIZE_DATAPOINTS
153)                 rndindex  = randperm(length(svmlabel));
154)                 svmdata   = svmdata(rndindex,:);
155)                 svmlabel  = svmlabel(rndindex);
156)             end

157)             

158)             decode = svm_single_crossval(svmlabel,svmdata,svmopts);
159)             % save the decode value to the corresponding coordinate

160) 
161)             coord= mapping.indexToCoordMap.get(iVoxel);

162)             x = coord(1);
163)             y = coord(2);
164)             z = coord(3);

165)             
166)             img3D(x,y,z) = ((decode/100)-0.5)*2; % range [-1:1]

167) 
168)         end  %for each voxel
169) 

170)         nii = make_nii(img3D,vdim,mask_image.hdr.hist.originator(1:3),16,sprintf('decode performance, time relative to onset: %g to %g sec',frameStart,frameEnd));
171)         save_nii(nii,fullfile(savePath,sprintf('%s-%03g',subjects{s}.name,timeShift)));

172)     
173)     end %for each timeslice

174)     display('rastering done');

175)     display(sprintf('result images saved to %s',savePath));
176)     

177)     elapsed{s}.decode = toc;

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

179) 
180)     %save memory!
181)     clear pst; 
182)     clear mask_image;

183) end
184) 
185) end