function output = runFBSImageMaskMode(header,subjectdata,fbsargs) addpath('NIFTI_20090325'); savePath = fullfile('output',datestr(now,30)); mkdir(savePath); args = header.args; subjects = subjectdata; nSubjects = numel(subjects); sessionlist = args.sessionList; pstopts = args.psthOpts; pstopts.eventList = args.eventList; pstopts.sessionList = sessionlist; radius = fbsargs.radius; timeline = header.timeline; timeline.frameShiftStart = header.frameShift.frameShiftStart; timeline.frameShiftEnd = header.frameShift.frameShiftEnd; timeline.decodeDuration = header.frameShift.decodeDuration; timePointArgs.labelMap = LabelMap(header.classDef.labelCells,header.classDef.conditionCells); timePointArgs.eventList = header.classDef.eventMatrix; timeLineStart = timeline.frameShiftStart; timeLineEnd = timeline.frameShiftEnd; RANDOMIZE_DATAPOINTS = header.svmrnd; svmopts = fbsargs.svmopts; disp(sprintf('batch processing %g subjects.',nSubjects)); for s = 1:nSubjects elapsed{s} = struct; %measure the timing disp(sprintf('processing subject %s.',subjects{s}.name)); % load image data disp('fetching volume definitions, please wait'); tic; volumes = spm_vol(getImageFileList(subjects{s}.dir,sessionlist,args.mask)); elapsed{s}.loadList = toc; disp('computing volume values, please wait'); tic [extr x y z] = calculateRoiImageData(volumes,subjects{s}.roiFile); elapsed{s}.calcExtr = toc; clear volumes; %save memory ?? disp('cleared volumes'); tic nVoxel = size(extr(1).dat,1); indexToCoordMap = java.util.HashMap; coordToIndexMap = java.util.HashMap; for iVoxel = 1:nVoxel coord = [x(iVoxel) y(iVoxel) z(iVoxel)] a = java.util.Vector(3,0); a.add(0,coord(1)); a.add(1,coord(2)); a.add(2,coord(3)); indexToCoordMap.put(iVoxel,coord); coordToIndexMap.put(a,iVoxel); end mapping.indexToCoordMap = indexToCoordMap; mapping.coordToIndexMap = coordToIndexMap; elapsed{s}.mapping = toc; % calculate psth warn = warning('off','all'); pstopts.des = subjects{s}.des; disp(sprintf('computing psth for %g voxel.',nVoxel)); tic for iVoxel = 1:nVoxel rawdata = []; for iImage = 1:length(extr); tmp = extr(iImage); rawdata = [rawdata tmp.dat(iVoxel)]; end % coordNameHelper = strcat('v',num2str(iVoxel)); pst{iVoxel} = calculatePST(args.timeline,pstopts,rawdata); % do not store in subjectStruct, so we can clear it later. % p = iVoxel/nVoxel*100; % if(mod(iVoxel,floor(nVoxel/100))==0) % sprintf(' %g%%\t complete',p); % end end elapsed{s}.psth = toc; disp('psth done'); warning(warn); clear extr; %run searchlight pause(0.001) % flush system event queue (respond to ctrl-c) tic display(sprintf('rastering %g coordinates with approx. %g mm radius',nVoxel,radius)); % for each timeslice globalStart = timeline.psthStart; globalEnd = timeline.psthEnd; decodeDuration = timeline.decodeDuration; labelMap = timePointArgs.labelMap; eventList = pstopts.eventList; res = timeline.tr*timeline.trFactor; % tmp = spm_imatrix(V(kImage).mat); % % vdim = tmp(7:9); % Voxel-Size mask_image = load_nii(subjects{s}.roiFile.fname); vdim = mask_image.hdr.dime.pixdim(2:4); if mask_image.hdr.dime.pixdim(1) == 1 vdim = vdim .* [-1 1 1]; end display(sprintf('starting timesliceing')); nSamplePoints = ((timeLineEnd-timeLineStart)/res) +1; if isempty( fbsargs.timeline ) fbsTimeLine = 1:nSamplePoints; fbsTimeLine = fbsTimeLine +globalStart; else fbsTimeLine = fbsargs.timeline; end for timeShiftIdx = fbsTimeLine % center timepoint && relative shift frameStartIdx = floor(-globalStart+1+timeShiftIdx - 0.5*decodeDuration); frameEndIdx = min(ceil(frameStartIdx+decodeDuration + 0.5*decodeDuration),-globalStart+globalEnd); img3D = zeros(size(mask_image)); %output image prepare for iVoxel = 1:nVoxel % linear structure avoids 3D-Loop. if (mod(iVoxel,100)== 0) display(sprintf('Status: %03u / %03u Timepoints, %05u / %05u Coordinates',find(fbsTimeLine == timeShiftIdx),length(fbsTimeLine),iVoxel,nVoxel)); pause(0.001) %flush system event queue end % get surrounding coordinate-IDs within radius sphere = fbs_buildSphere(mapping,iVoxel,radius,vdim); % pause %build svm inputmatrix svmdata = []; svmlabel = []; anyvoxel = 1; for pstConditionGroup = 1:size(pst{1,anyvoxel},2) for dp = 1:size(pst{1,anyvoxel}{1,pstConditionGroup},1) % data point svmlabel = [svmlabel; lm_getSVMLabel(labelMap,eventList(pstConditionGroup,1))]; row = []; for voxel = 1:size(sphere,2) row = [row, pst{1,voxel}{1,pstConditionGroup}(dp,frameStartIdx:frameEndIdx)]; % label,[value,value,...],[value,value,...]... end svmdata = [svmdata; row]; end end if RANDOMIZE_DATAPOINTS % [svmlabel svmdata] rndindex = randperm(length(svmlabel)); svmdata = svmdata(rndindex,:); svmlabel = svmlabel(rndindex); end decode = svm_single_crossval(svmlabel,svmdata,svmopts); % save the decode value to the corresponding coordinate coord = mapping.indexToCoordMap.get(iVoxel) x = coord(1); y = coord(2); z = coord(3); img3D(x,y,z) = ((decode/100)-0.5)*2; % range [-1:1] end %for each voxel nii = make_nii(img3D,vdim,mask_image.hdr.hist.originator(1:3),16,... sprintf('decode performance, time relative to onset: %g to %g sec',frameStartIdx,frameEndIdx)); save_nii(nii,fullfile(savePath,sprintf('%s-%03g',subjects{s}.name,timeShiftIdx))); end %for each timeslice display('rastering done'); display(sprintf('result images saved to %s',savePath)); elapsed{s}.decode = toc; assignin('base','timing',elapsed); %save memory! clear pst; clear mask_image; end end