442 lines
18 KiB
Matlab
442 lines
18 KiB
Matlab
classdef BSM<handle
|
|
%BSM Summary of this class goes here
|
|
% Detailed explanation goes here
|
|
|
|
properties %public properties
|
|
SPMs={}; %field to contain the full path of all SPM.mat files to be analysed
|
|
Acfl %columns in SPM.xX.X to ignore
|
|
Type=4; %BSM type 1..4 5: none %default is option 4
|
|
KeepMP = false; %retain Motion paramters
|
|
Prefix = 'bs'; %beta series prefix
|
|
Clean = true;
|
|
end
|
|
properties (Access = private)%private properties
|
|
SPM_orig %temporary storage for SPM.mat file
|
|
ActiveColumns %store the columns to be analized
|
|
U_orig %store the U field of SPM
|
|
TargetFolder %target folder
|
|
% option below obsolete?
|
|
BSMToolboxPath%store path to BSMToolbox
|
|
SPMPath;%store path to SPM folder
|
|
SPMIndexList=[];
|
|
MonitorState=false;
|
|
end
|
|
properties (Dependent) %dependent properties
|
|
NumSPMs% number of SPM.mat files to process
|
|
|
|
end
|
|
properties (Constant = true)%constants
|
|
BSMMethod={'Traditional BSM Rissman et al', 'Adapted BSM Rissman et al', 'Fast approach (Mumford et al)', 'Adapted fast approach','none'};
|
|
end
|
|
|
|
methods %user callable methods
|
|
function Monitor(obj,varargin)
|
|
if numel(varargin)>1
|
|
id = 'BSM:Monitor:TooManyInputs';
|
|
msg= sprintf('Too many inputs given.\n Use: \t obj.Monitor\n\tobj.Monitor(true)\n\tobj.monitor(false)\n');
|
|
throw(MException(id,msg))
|
|
end
|
|
if numel(varargin)==0
|
|
obj.MonitorState=xor(obj.MonitorState,true);
|
|
else
|
|
if ~islogical(varargin{1})
|
|
id = 'BSM:Monitor:InputIsNotLogical';
|
|
msg= sprintf('Input is not of type logical\n');
|
|
throw(MException(id,msg))
|
|
end
|
|
obj.MonitorState=varargin{1};
|
|
end
|
|
end
|
|
function Calc(obj,varargin)
|
|
% obj.calc to apply calculate to all SPM.mat files
|
|
% obj.calc(n) to calculated the nth SPM.mat file
|
|
% obj.calc([1 3 6]) to indicate which SPM.mat file to calculate
|
|
|
|
|
|
SPMidx=[];% will contain the index to the SPMs to work on
|
|
if logical(numel(varargin))
|
|
SPMidx=varargin{1};
|
|
end
|
|
%check number of input arguments
|
|
if ~logical(obj.NumSPMs)%nothing to do
|
|
fprintf('== == == Nothing to do: done == == ==\n')
|
|
return;
|
|
end
|
|
if isempty(SPMidx)
|
|
SPMidx=1:obj.NumSPMs;
|
|
end
|
|
if ~isnumeric(SPMidx)||~all(isfinite(SPMidx))||min(SPMidx)<1||max(SPMidx)>obj.NumSPMs
|
|
id='BSM:Calc:InvalidIndex';
|
|
msg='Invalid Index';
|
|
throw(MException(id,msg));
|
|
end
|
|
if numel(SPMidx)>1
|
|
obj.SPMIndexList=SPMidx;
|
|
SPMidx=obj.SPMIndexList(end);
|
|
try
|
|
obj.Calc(SPMidx); %handle the single case
|
|
catch ME
|
|
%if the identifier matches, just continue with next SPM
|
|
if ~any(strcmpi(ME.identifier,...
|
|
{'BSM:CoreCalc:SPMxBFOrder:invalid'...
|
|
,'BSM:CoreCalc:SPMSess:invalid'...
|
|
}))
|
|
rethrow(ME)
|
|
else
|
|
disp(ME.message);
|
|
end
|
|
end
|
|
obj.SPMIndexList=obj.SPMIndexList(1:end-1);
|
|
if isempty(obj.SPMIndexList)
|
|
return;
|
|
else
|
|
obj.Calc(obj.SPMIndexList);% recursively run the rest of the list
|
|
end
|
|
else % deal with single entry
|
|
obj.CoreCalc(SPMidx)
|
|
end
|
|
|
|
|
|
end
|
|
end
|
|
methods %constructor and set/get methods
|
|
function obj = BSM(varargin)
|
|
%BSM Construct an instance of this class
|
|
%option: parse the info of a matlab batch file? see
|
|
%tbx_cf_bsm?
|
|
%TODO:
|
|
% set default values
|
|
obj.SPMPath=[]; %get the SPM path
|
|
if nargin>0
|
|
[val,TF]=obj.MyGetVal('SPMs',varargin);if TF;obj.SPMs=val;end
|
|
[val,TF]=obj.MyGetVal('Type',varargin);if TF;obj.Type=val;end
|
|
[val,TF]=obj.MyGetVal('Acfl',varargin);if TF;obj.Acfl=val;end
|
|
end
|
|
% initialize spm job enviourment
|
|
spm_jobman('initcfg');
|
|
end
|
|
function val=get.NumSPMs(obj)
|
|
val=numel(obj.SPMs);
|
|
end
|
|
function set.SPMPath(obj,~)
|
|
SPMFound=which('spm'); %check if SPM in searchpath
|
|
if isempty(SPMFound) %SPM not found
|
|
error('please add SPM to your path before trying to run BSM');
|
|
%to do add a check on spm version
|
|
end
|
|
obj.SPMPath=fileparts(SPMFound);
|
|
end
|
|
function set.Clean(obj,varargin)
|
|
if isempty(varargin)
|
|
%% this option is obsolete
|
|
obj.Clean=xor(obj.Clean,true); %togle state
|
|
elseif ~islogical(varargin{1})
|
|
id = 'BSM:SetClean:InvalidType';
|
|
msg= 'Value passed to obj.Clean is not of type logical';
|
|
throw(MException(id,msg));
|
|
else
|
|
obj.Clean=varargin{1};
|
|
end
|
|
end
|
|
function set.Prefix(obj,val)
|
|
if ischar(val)
|
|
obj.Prefix=val;
|
|
else
|
|
id='BSM:SetPrefix:InvalidInput';
|
|
msg='Input needs to be a char array';
|
|
throw(MException(id,msg));
|
|
end
|
|
end
|
|
function set.Type(obj,val)
|
|
%check allowed input
|
|
if ~isnumeric(val)|| any(~ismember(val,1:5))
|
|
id='BSM:Type:invalid';
|
|
msg='use a value in the range 1..5';
|
|
throw(MException(id,msg));
|
|
end
|
|
obj.Type=val;
|
|
end
|
|
function set.Acfl(obj,val)
|
|
%allow clearing
|
|
if isempty(val)
|
|
obj.Acfl=[];
|
|
end
|
|
if ~isnumeric(val) || any(val<1) || any(~isfinite(val))
|
|
id='BSM:Acfl:invalid';
|
|
msg='set an array of positive integer values to the Acfl field';
|
|
throw(MException(id,msg));
|
|
end
|
|
obj.Acfl=val;
|
|
|
|
end
|
|
function set.SPMs(obj,val)
|
|
%set the SPMs to be analysed
|
|
%obj.SPMs=<val>
|
|
%possible values for <val>
|
|
% {} : clear the SPMs field
|
|
% {'ui'} : ask user input
|
|
if ~iscell(val)
|
|
fprintf('%s\n',...
|
|
'Use any of the following options',...
|
|
'obj.SPMs={}',...
|
|
'obj.SPMs={''ui''}',...
|
|
'obj.SPMs={<Ful path to SPM.mat>}');
|
|
return
|
|
end
|
|
if isempty(val) %allow empty input
|
|
return;
|
|
end
|
|
if strcmpi(val{1},'ui')
|
|
[P,sts]=spm_select([1 Inf],'mat','select the SPM.mat files',obj.SPMs,pwd,'^SPM.mat');
|
|
if ~logical(sts);return;end
|
|
if iscell(P)
|
|
val=P;
|
|
else
|
|
val=cellstr(P);
|
|
end
|
|
end
|
|
% capture file not found somewhat clever
|
|
TF=cellfun(@exist,val);
|
|
if any(not(TF))
|
|
id='BSM:SetSPMs:FileNotFound';
|
|
msg=sprintf('file not found\t%s\n',val{not(TF)});
|
|
throw(MException(id,msg));
|
|
end
|
|
obj.SPMs=val;
|
|
end
|
|
end
|
|
methods (Access = protected)%private functions
|
|
function CoreCalc(obj,idx)
|
|
fprintf('working on :%d: %s\n',idx,obj.SPMs{idx});
|
|
% how to store the current status of warning and restore curent
|
|
WarnStat=warning;
|
|
warning('OFF', 'MATLAB:MKDIR:DirectoryExists');
|
|
%this try catch block will handle any errors that might ocure
|
|
%warning will not remain off by accident
|
|
try
|
|
SpmFile=obj.SPMs{idx};
|
|
load(SpmFile); %#ok<LOAD>
|
|
%% exception block
|
|
%NB these exception has not been tested properly by the
|
|
if (SPM.xBF.order ~=1 ) %#ok<*NODEF>
|
|
%TestBSMClass program
|
|
%throw an error. The calling function should properly
|
|
%handle the error.
|
|
id='BSM:CoreCalc:SPMxBFOrder:invalid';
|
|
msg=sprintf('%s\n','Problem found', ...
|
|
'No canonical HRF was used to model the data', ...
|
|
'Current file will be skipped');
|
|
throw(MException(id,msg));
|
|
end
|
|
if (length(SPM.Sess) > 1 )
|
|
id='BSM:CoreCalc:SPMSess:invalid';
|
|
msg=sprintf('%s\n',...
|
|
'Problem found', ...
|
|
'Multiple sessions in one SPM.mat not supported', ...
|
|
'Current file will be skipped');
|
|
throw(MException(id,msg));
|
|
end
|
|
% end exception block
|
|
%% set local variables
|
|
SPMOrig = SPM;
|
|
UOrig = SPM.Sess(1).U;
|
|
[ActCol,~,~]=obj.GetColDef(UOrig);
|
|
LocalTargetFolder = fullfile(fileparts(SpmFile),'BSM');%RR replace BSM by a obj.BSMDir so user can modify this
|
|
%% create directory
|
|
mkdir(LocalTargetFolder);
|
|
% loop over ActCol
|
|
for c=1:length(ActCol)
|
|
% get current col, get other col
|
|
CurCol=ActCol(c);
|
|
OthCol=setdiff(1:length(ActCol),CurCol);
|
|
TargetFolderCurCol=fullfile(LocalTargetFolder,sprintf('col_%04d', CurCol));
|
|
NEvents=numel(UOrig(CurCol).ons);
|
|
DirTrialList=cell(NEvents,1);
|
|
BetaTrialList=cell(NEvents,1);
|
|
for e=1:NEvents% loop over events
|
|
U=obj.GetU(UOrig,e,CurCol,OthCol);%build U
|
|
DirTrialList{e}=obj.GenSPM(U,SPMOrig,TargetFolderCurCol,e);%build SPM
|
|
obj.RunJob(DirTrialList{e});%process SPM
|
|
BetaTrialList{e}=spm_select('ExtFPList',DirTrialList{e},'^beta_0001.*.nii$',1);
|
|
end% end loop
|
|
obj.PackBetas(BetaTrialList,LocalTargetFolder,CurCol)%pack per event result
|
|
obj.CleanUp(TargetFolderCurCol);
|
|
end % end loop act. col
|
|
% restore warning status
|
|
warning(WarnStat);
|
|
catch ME %turn matlab warning back on
|
|
warning(WarnStat);
|
|
rethrow(ME);
|
|
end
|
|
end
|
|
function CleanUp(~,D)
|
|
if obj.Clean
|
|
try
|
|
rmdir(D,'s');
|
|
catch ME
|
|
disp(ME)
|
|
id = 'BSM:Clean:Failed';
|
|
msg= 'Cleaning failed';
|
|
throw(MException(id,msg));
|
|
end
|
|
end
|
|
end
|
|
function PackBetas(obj,P,wd,Col)
|
|
FPFout=fullfile(wd,sprintf('%s_%04d.nii',obj.Prefix,Col));
|
|
matlabbatch{1}.spm.util.cat.vols = P;
|
|
matlabbatch{1}.spm.util.cat.name = FPFout;
|
|
matlabbatch{1}.spm.util.cat.dtype = 4;
|
|
spm_jobman('run', matlabbatch);
|
|
end
|
|
|
|
function RunJob(~,swd)
|
|
%% estimate model:
|
|
matlabbatch{1}.spm.stats.fmri_est.spmmat = {fullfile(swd,'SPM.mat')};
|
|
matlabbatch{1}.spm.stats.fmri_est.method.Classical = 1;
|
|
spm_jobman('run', matlabbatch);
|
|
end
|
|
function swd=GenSPM(obj,U,SPMOrig,TargetFolderCol,Event)
|
|
SPM=SPMOrig;
|
|
TargetFolderColTrial=fullfile(TargetFolderCol,sprintf('trial_%04d',Event));
|
|
mkdir(TargetFolderColTrial);
|
|
SPM.swd = TargetFolderColTrial;
|
|
SPM.Sess.U=U;
|
|
U=spm_get_ons(SPM,1); % to have SPM autocomplete the U structure
|
|
[X,~,Fc]=spm_Volterra(U,SPM.xBF.bf,SPM.xBF.Volterra);
|
|
|
|
k=SPM.nscan(1);
|
|
fMRI_T = SPM.xBF.T;
|
|
fMRI_T0 = SPM.xBF.T0;
|
|
|
|
X = X((0:(k-1))*fMRI_T +fMRI_T0 + 32,:);%@JB what does the +32 do?
|
|
%orthogonalize
|
|
for k=1:numel(Fc)
|
|
X(:,Fc(k).i)=spm_orth(X(:,Fc(k).i));
|
|
end
|
|
|
|
Xname=cell(1,size(X,2));
|
|
for k = 1:size(X,2)
|
|
Xname{k} = [sprintf('Sn(%i) ',1) SPM.Sess(1).U(k).name{1}];
|
|
end
|
|
|
|
%% reinstate motion parameters
|
|
if (obj.KeepMP)&& ~isempty(SPMOrig.Sess(1).C.C)
|
|
C = SPMOrig.Sess(1).C.C;
|
|
Cname = job.SPM_orig.Sess(1).C.name;
|
|
|
|
X = cat(2,X,spm_detrend(C));
|
|
Xname = cat(2,Xname, Cname(:));
|
|
end
|
|
%% reinstate baseline column
|
|
X(:,end+1) = 1;
|
|
Xname = cat(2,Xname,{sprintf('Sn(%i) %s',1,'constant')});
|
|
%% monitor
|
|
if obj.MonitorState
|
|
figure(100);
|
|
clf;
|
|
imagesc(X);
|
|
pause(0.5);
|
|
end
|
|
%% Rebuild SPM
|
|
SPM.xX.X = X;
|
|
SPM.xX.iB = size(X,2);
|
|
SPM.xX.iC = 1:size(X,2)-1;
|
|
SPM.xX.name = Xname;
|
|
try
|
|
SPM.xX = rmfield(SPM.xX,'W');%remove prewhitening matrix
|
|
catch
|
|
end
|
|
save(fullfile(SPM.swd,'SPM.mat'),'SPM');
|
|
swd=SPM.swd;
|
|
end
|
|
function U=GetU(obj,UOrig,Event,Col,OCol)
|
|
U=struct('name',{},'ons', [], 'dur', [],'orth', [], 'P', [], 'dt', [], 'u', [],'pst',[]);
|
|
idx=1;
|
|
%split out current event as first condition
|
|
U(idx).name{1} = sprintf('%s_%d',UOrig(Col).name{1},Event);
|
|
U(idx).ons = UOrig(Col).ons(Event);
|
|
U(idx).dur = UOrig(Col).dur(Event);
|
|
U(idx).orth = UOrig(Col).orth;
|
|
U(idx).P = UOrig(Col).P;
|
|
U(idx).dt = UOrig(Col).dt;
|
|
% u
|
|
% pst
|
|
%optional add all other events as the next condition
|
|
switch obj.Type
|
|
case {1,2}
|
|
idx = idx +1;
|
|
OtherTrials=setdiff(1:numel(UOrig(Col).ons),Event);
|
|
U(idx).name{1} = sprintf('%s_oe',UOrig(Col).name{1});
|
|
U(idx).ons = UOrig(Col).ons(OtherTrials);
|
|
U(idx).dur = UOrig(Col).dur(OtherTrials);
|
|
U(idx).orth = UOrig(Col).orth;
|
|
U(idx).P = UOrig(Col).P;
|
|
U(idx).dt = UOrig(Col).dt;
|
|
% u
|
|
% pst
|
|
case {3,4}
|
|
otherwise
|
|
end
|
|
% deal with other conditions
|
|
if ~isempty(OCol)
|
|
switch obj.Type
|
|
case {1,4} %put all other conditions in one column
|
|
idx = idx + 1;
|
|
U(idx).name{1} = 'AllOtherConditions';
|
|
U(idx).ons = cat(1,UOrig(OCol).ons);
|
|
for k = OCol
|
|
if numel(UOrig(k).dur)==1
|
|
dur=repmat(UOrig(k).dur,size(UOrig(k).ons));
|
|
elseif numel(UOrig(k).dur)==numel(UOrig(k).ons)
|
|
dur=UOrig(k).dur;
|
|
else
|
|
id='BSM:GetU:InvalidNumberDuration';
|
|
msg=sprintf('%s\n'...
|
|
,'Number of durations is neither one nor the number of onsets'...
|
|
);
|
|
throw(MException(id,msg));
|
|
end
|
|
U(idx).dur=cat(1,U(idx).dur,dur);
|
|
U(idx).orth = UOrig(k).orth;
|
|
U(idx).P = UOrig(k).P;
|
|
U(idx).dt = UOrig(k).dt;
|
|
end
|
|
case {2,3} %keep other conditions each in their own column
|
|
U=cat(2,U,UOrig(OCol));
|
|
otherwise
|
|
end
|
|
%clear out fields u and pst They will be rebuild when SPM
|
|
%is created
|
|
[U(:).pst]=deal([]);
|
|
[U(:).u]=deal([]);
|
|
end
|
|
end
|
|
function [ActCol,PasCol,AllCol]=GetColDef(obj,U)
|
|
AllCol=1:length(U);
|
|
if ~isempty(obj.Acfl)
|
|
PasCol = AllCol(obj.Acfl); %cf line 53 in spm_bsm.m
|
|
else
|
|
PasCol = [];
|
|
end
|
|
ActCol = setdiff(AllCol,PasCol);
|
|
end
|
|
function[val,TF]=MyGetVal(~,selector,List)
|
|
TF=strcmpi(selector,List);
|
|
if any(TF) %spmpath is set as option
|
|
k=find(TF, 1 );%find the first instance
|
|
if (k+1)>numel(List)
|
|
id='BSM:getval:ReadBeyondEndInput';
|
|
msg=sprintf('use name value pairs\nlast entry is a name');
|
|
throw(MException(id,msg));
|
|
end
|
|
val=List{k+1};
|
|
TF=true;
|
|
else
|
|
val=false;
|
|
TF=false;
|
|
end
|
|
end
|
|
end
|
|
end
|
|
|