Commit 2cbd4786 authored by Stephen Smith's avatar Stephen Smith Committed by Duncan Mortimer
Browse files

Updates to fix_3_clean.m

parent f284a66e
%
% fix_3_clean(fixlist,aggressive,domot,hp) - apply the FIX cleanup to filtered_func_data
% fix_3_clean(fixlist,aggressive,domot,hp,0) - CIFTI processing only (do not process volumetric data)
%
% fixlist is a vector of which ICA components to remove (starting at 1 not 0)
%
......@@ -13,7 +14,7 @@
% hp>0 the fullwidth (2*sigma) of the highpass, in seconds (not TRs)
%
function fix_3_clean(fixlist,aggressive,domot,hp)
function fix_3_clean(fixlist,aggressive,domot,hp,varargin)
if (isdeployed)
aggressive = str2num(aggressive)
domot = str2num(domot)
......@@ -26,6 +27,11 @@ WBC=getenv('FSL_FIX_WBC');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
DOvol=1;
if length(varargin) > 0 && varargin{1} == 0
DOvol=0;
end
%%%% read set of bad components
DDremove=load(fixlist);
......@@ -39,21 +45,25 @@ if exist('Atlas.dtseries.nii','file') == 2
path(path,CIFTI);
BO=ciftiopen('Atlas.dtseries.nii',WBC);
if hp==0
BO.cdata=detrend(BO.cdata')';
meanBO=mean(BO.cdata,2);
BO.cdata=detrend(BO.cdata')'; BO.cdata=BO.cdata+repmat(meanBO,1,size(BO.cdata,2));
end
if hp>0
BOdimX=size(BO.cdata,1); BOdimZnew=ceil(BOdimX/100); BOdimT=size(BO.cdata,2);
meanBO=mean(BO.cdata,2); BO.cdata=BO.cdata-repmat(meanBO,1,size(BO.cdata,2));
save_avw(reshape([BO.cdata ; zeros(100*BOdimZnew-BOdimX,BOdimT)],10,10,BOdimZnew,BOdimT),'Atlas','f',[1 1 1 TR]);
call_fsl(sprintf('fslmaths Atlas -bptf %f -1 Atlas',0.5*hp/TR));
grot=reshape(read_avw('Atlas'),100*BOdimZnew,BOdimT); BO.cdata=grot(1:BOdimX,:); clear grot;
grot=reshape(read_avw('Atlas'),100*BOdimZnew,BOdimT); BO.cdata=grot(1:BOdimX,:); clear grot; BO.cdata=BO.cdata+repmat(meanBO,1,size(BO.cdata,2));
ciftisave(BO,'Atlas_hp_preclean.dtseries.nii',WBC); % save out noncleaned hp-filtered data for future reference, as brainordinates file
end
end
%%%% read NIFTI version of the data
cts=read_avw('filtered_func_data');
ctsX=size(cts,1); ctsY=size(cts,2); ctsZ=size(cts,3); ctsT=size(cts,4);
cts=reshape(cts,ctsX*ctsY*ctsZ,ctsT)';
if DOvol
cts=read_avw('filtered_func_data');
ctsX=size(cts,1); ctsY=size(cts,2); ctsZ=size(cts,3); ctsT=size(cts,4);
cts=reshape(cts,ctsX*ctsY*ctsZ,ctsT)';
end
%%%% read and prepare motion confounds
confounds=[];
......@@ -68,8 +78,10 @@ ICA=functionnormalise(load(sprintf('filtered_func_data.ica/melodic_mix')));
if aggressive == 1
sprintf('aggressive cleanup')
confounds=[confounds ICA(:,DDremove)];
cts = cts - (confounds * (pinv(confounds,1e-6) * cts));
if DObrainord == 1
if DOvol
cts = cts - (confounds * (pinv(confounds,1e-6) * cts));
end
if DObrainord
BO.cdata = BO.cdata - (confounds * (pinv(confounds,1e-6) * BO.cdata'))';
end
else
......@@ -77,23 +89,45 @@ else
if domot == 1
% aggressively regress out motion parameters from ICA and from data
ICA = ICA - (confounds * (pinv(confounds,1e-6) * ICA));
cts = cts - (confounds * (pinv(confounds,1e-6) * cts));
if DObrainord == 1
if DOvol
cts = cts - (confounds * (pinv(confounds,1e-6) * cts));
end
if DObrainord
BO.cdata = BO.cdata - (confounds * (pinv(confounds,1e-6) * BO.cdata'))';
end
end
betaICA = pinv(ICA,1e-6) * cts; % beta for ICA (good *and* bad)
cts = cts - (ICA(:,DDremove) * betaICA(DDremove,:)); % cleanup
if DObrainord == 1
betaICA = pinv(ICA,1e-6) * BO.cdata'; % beta for ICA (good *and* bad)
if DOvol
betaICA = pinv(ICA,1e-6) * cts; % beta for ICA (good *and* bad)
cts = cts - (ICA(:,DDremove) * betaICA(DDremove,:)); % cleanup
end
if DObrainord
betaICA = pinv(ICA,1e-6) * BO.cdata'; % beta for ICA (good *and* bad)
BO.cdata = BO.cdata - (ICA(:,DDremove) * betaICA(DDremove,:))'; % cleanup
end
end
%%%% save cleaned data to file
save_avw(reshape(cts',ctsX,ctsY,ctsZ,ctsT),'filtered_func_data_clean','f',[1 1 1 1]);
call_fsl('fslcpgeom filtered_func_data filtered_func_data_clean');
if DObrainord == 1
if DOvol
save_avw(reshape(cts',ctsX,ctsY,ctsZ,ctsT),'filtered_func_data_clean','f',[1 1 1 1]);
call_fsl('fslcpgeom filtered_func_data filtered_func_data_clean');
end
if DObrainord
ciftisave(BO,'Atlas_clean.dtseries.nii',WBC);
end
%%%% compute variance normalization field and save to file
DDSignal=setdiff([1:1:size(ICA,2)],DDremove); % select signal components
if DOvol
betaICA = pinv(ICA(:,DDSignal),1e-6) * cts; % beta for ICA (good, bad already removed)
cts = cts - (ICA(:,DDSignal) * betaICA); % remove signal components to compute unstructured noise timeseries
cts = std(cts); % compute variance normalization map
save_avw(reshape(cts',ctsX,ctsY,ctsZ),'filtered_func_data_clean_vn','f',[1 1 1 1]);
call_fsl('fslcpgeom filtered_func_data filtered_func_data_clean_vn -d');
end
if DObrainord
betaICA = pinv(ICA(:,DDSignal),1e-6) * BO.cdata'; % beta for ICA (good, bad already removed)
BO.cdata = BO.cdata - (ICA(:,DDSignal) * betaICA)'; % remove signal components to compute unstructured noise timeseries
BO.cdata = std(BO.cdata,[],2); % compute variance normalization map
ciftisavereset(BO,'Atlas_clean_vn.dscalar.nii',WBC);
end
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment