private/runFBSImageMaskMode.m (faae5fb) - SVMCrossVal.git

runFBSImageMaskMode.m

function output =  runFBSImageMaskMode(header,subjectdata,fbsargs)
addpath(fullfile(getTbxPath,'NIFTI_20090325'));

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

subjectResultDir = subjectdata{1}.dir;
savePath = fullfile(subjectResultDir,'results','fbs',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;

radiusList = 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
    if isDebug
        figure;
        hold on;
        for i = 1:size(pst,2)
            plot(mean(pst{:,i}{1}),'r-');
            plot(mean(pst{:,i}{2}),'b-');
        end
        hold off;
    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',nVoxel));
    % 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);
    tmp = strfind(mask_image.fileprefix,filesep);
    maskname = mask_image.fileprefix(tmp(end)+1:end);
    
    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 
    %preferred!
	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);

for rIdx = 1:length(radiusList)
            img3D{rIdx}(:,:,:) = zeros(size(mask_image)); %output image prepare
        end
        
        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
            for rIdx  = 1: length(radiusList)
                radius = radiusList(rIdx);
                % get surrounding coordinate-IDs within radius
                sphere = fbs_buildSphere(mapping,iVoxel,radius,vdim);
%                 for i = 1: length(sphere)
%                      sphere(i )
%                      mapping.indexToCoordMap.get(sphere(i))
%                 end
                
                %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 = sphere
                            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);
                
                if NORM_DECODE
                    img3D{rIdx}(x,y,z) = ((decode/100)-0.5)*2; % range [-1:1]
                else
                    img3D{rIdx}(x,y,z) = (decode/100); % range [0:1]
                end

end %for each radius
        end %for each voxel

for rIdx = 1:length(radiusList)
            radius = radiusList(rIdx);
            nii = make_nii(img3D{rIdx},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-%s-r%g-t%+03g',subjects{s}.name,maskname,radius,timeShiftIdx)));
        end
    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

SVMCrossVal.git

fbs save path, axels stimolos.mat