Homec4science
Diffusion Slepians (master) Browse

Slepi/Code/Utilities/MyPreprocessing/ntimeCourseExtractor_TB.m
22f22bf83060master

Last Change
Raw File

ntimeCourseExtractor_TB.m
View Options

%% Generates atlased time courses for five atlas parcellations: AAL, Craddock
% (300 and 900 ROIs), Shen and James atlases
function [TCS_A1,TCS_A2,TCS_A3,TCS_A4,TCS_A5]=ntimeCourseExtractor_TB(paths,pre_f,pre_f2,reg,Ai,A,codeBook,scansIdx,fid,varargin)
%% 1. Default options
[DEBUGMODE,blockDriftCorrection,firstEV,regressOutCSF,regressOutWM,regressOutCO2,regressOutMot,Mot12,...
removeFirstNscans,doBootstrap,boot_NN,boot_nReps,tryParsingFileOrder,...
Detrend,Band,TR,Despike,ScrubbingMethod]=process_options_TB(varargin,...
'DEBUGMODE',false,'blockDriftCorrection','none','firstEV',false,...
'regressOutCSF',true,'regressOutWM',1,'regressOutCO2',1,'regressOutMot',true,'Mot12',0,'removeFirstNscans',10,'doBootstrap',false,...
'boot_NN',10,'boot_nReps',100,'tryParsingFileOrder',true,'Detrend',1,...
'Band',[0.01 0.15],'TR',[],'Despike',1,'ScrubbingMethod',[]);
%% 2. Information on one functional volume
% Gets one volume only
[fMeanFN_noPath, foosubdirs] = spm_select('List',paths.f,['^' pre_f '.*\.nii$']);
fMeanFN_noPath=fMeanFN_noPath(1+removeFirstNscans,:);
% Gets the path and the volume itself
fMean_fn=fullfile(paths.f,fMeanFN_noPath);
rVol=spm_read_vols(spm_vol(fMean_fn));
%% 3. LOAD CSF MASK
if regressOutCSF==true
fprintf('Extracting CSF...');
if isempty(strfind(pre_f,reg))
[segFN_noPath, foosubdirs] = spm_select('List',paths.sseg,'^wc3.*\.nii$'); % new file with overlap with mask - in individual space case
sSeg_fn=fullfile(paths.sseg,segFN_noPath);
CSFV=mapVolumeToVolume_TB(sSeg_fn,fMean_fn);
CSFmaskLidx=CSFV;
clear CSFV;
end
fprintf('\n');
end
%% 4. LOAD WM MASK
if regressOutWM==true
fprintf('Extracting WM...')
if isempty(strfind(pre_f,reg)) %individual space
[segFN_noPath, foosubdirs] = spm_select('List',paths.sseg,'^wc2.*\.nii$');
sSeg_fn=fullfile(paths.sseg,segFN_noPath);
WMV=mapVolumeToVolume_TB(sSeg_fn,fMean_fn);
WMmaskLidx=WMV;
clear WMV;
end
fprintf('\n');
end
%% 5. EXTRACT FMRI TIMECOURSES
disp('Extracting time courses...');
% Names of all the functional files
fVolsFNlist=getImageFNinAcqOrder_TB(paths.f,pre_f,'nii',[]);
fVolsFNlist_VX=getImageFNinAcqOrder_TB(paths.f,pre_f2,'nii',[]);
fVolsFNlist_fullpath=cellfun(@(x) fullfile(paths.f,x), fVolsFNlist,'UniformOutput',false);
fVolsFNlist_fullpath_VX=cellfun(@(x) fullfile(paths.f,x), fVolsFNlist_VX,'UniformOutput',false);
% read all headers and files in temporal order into a 4-D array
V0i=spm_vol(fVolsFNlist_fullpath);
V0i_VX=spm_vol(fVolsFNlist_fullpath_VX);
if removeFirstNscans~=0
V0i=V0i((removeFirstNscans+1):end);
V0i_VX=V0i_VX((removeFirstNscans+1):end);
end
% Indices of the frames of interest (common to smoothed and unsmoothed
% cases
V0idx=1:length(V0i);
% Creates the 4D volumes
V0=zeros(V0i{1}.dim(1),V0i{1}.dim(2),V0i{1}.dim(3),length(V0i),'single');
V0_VX=zeros(V0i{1}.dim(1),V0i{1}.dim(2),V0i{1}.dim(3),length(V0i),'single');
for iter=1:length(V0i)
V0(:,:,:,V0idx(iter))=spm_read_vols(V0i{iter});
V0_VX(:,:,:,V0idx(iter))=spm_read_vols(V0i_VX{iter});
end
%% 6. Detrending
disp('Detrending...');
AllVolume=reshape(V0,[],size(V0,4))'; % time x space
AllVolume_VX=reshape(V0_VX,[],size(V0_VX,4))'; % time x space
CUTNUMBER=10; % cut space in parts
if Detrend==1
SegmentLength = ceil(size(AllVolume,2) / CUTNUMBER);
for iCut=1:CUTNUMBER
if iCut~=CUTNUMBER
Segment = (iCut-1)*SegmentLength+1 : iCut*SegmentLength;
else
Segment = (iCut-1)*SegmentLength+1 : size(AllVolume,2);
end
if Detrend
AllVolume(:,Segment) = detrend(AllVolume(:,Segment));
AllVolume_VX(:,Segment) = detrend(AllVolume_VX(:,Segment));
end
end
end
V0=reshape(AllVolume',size(V0,1),size(V0,2),size(V0,3),size(V0,4));
V0_VX=reshape(AllVolume_VX',size(V0_VX,1),size(V0_VX,2),size(V0_VX,3),size(V0_VX,4));
% Sets NaN elements of rVol to 0
rVol(isnan(rVol)) = 0;
%% 7. COMPUTE AVG CSF TIMECOURSE
if regressOutCSF==true
disp(' Computing average CSF signal...');
% Sets to 0 the elements that are null in the functional images
CSFmaskLidx(~logical(rVol))=0;
% Returns the locations in X, Y and Z of the elements from the mask
[Ci,Cj,Ck]=ind2sub(size(V0(:,:,:,1)),find(CSFmaskLidx));
% How many voxels
nCSFvoxels=numel(Ci);
WriteInformation_TB(fid,['Voxels CSF: ',num2str(nCSFvoxels),'\t']);
% Creates the time course itself
tcCSF=zeros(1,1,1,size(V0,4));
for idx=1:numel(Ci)
tcCSF=tcCSF+V0(Ci(idx),Cj(idx),Ck(idx),:);
end
CSFavg=squeeze(tcCSF/nCSFvoxels)';
CSFavg_demean=CSFavg-mean(CSFavg);
clear tcCSF;
end
%% 8. COMPUTE AVG WM TIMECOURSE
if regressOutWM==true
disp(' Computing average WM signal...');
WMmaskLidx(~logical(rVol))=0;
[Ci,Cj,Ck]=ind2sub(size(V0(:,:,:,1)),find(WMmaskLidx));
nWMvoxels=numel(Ci); % also, sum(CSFmaskLidx)
WriteInformation_TB(fid,['WM: ',num2str(nWMvoxels),'\t']);
tcWM=zeros(1,1,1,size(V0,4));
for idx=1:numel(Ci)
tcWM=tcWM+V0(Ci(idx),Cj(idx),Ck(idx),:);
end
WMavg=squeeze(tcWM/nWMvoxels)';
WMavg_demean=WMavg-mean(WMavg);
clear tcWM;
end
%% 9. Construction of the matrix of regressors
fprintf('Regressing out nuisance covariates...');
% Constant, linear and quadratic regressors
X=[ones(numel(CSFavg),1) [1:numel(CSFavg)]'/numel(CSFavg) ...
[1:numel(CSFavg)].^2'/(numel(CSFavg)^2)];
if regressOutWM
% Adds the WM regressor
X=[X,WMavg_demean'];
end
if regressOutCSF
% Adds the CSF regressor
X=[X,CSFavg_demean'];
end
% If we want to regress motion, we look for the text file and sample
% the motion columns
if regressOutMot
mot = dir(fullfile(paths.f,'rp*.txt'));
Cov=load(fullfile(paths.f,mot(1).name));
X=[X,Cov(removeFirstNscans+1:end,1:6)]; %take first 6 (not derivative if they are there)
end
%% 10. Regression itself
disp('Regressing out covariates...');
% Regression for both the data to atlas and the voxelwise data
V0Cov = zeros(size(V0));
V0Cov_VX =zeros(size(V0_VX));
for i=1:size(V0,1)
fprintf('.');
for j=1:size(V0,2)
for k=1:size(V0,3)
% We keep the residual from the regression process (i.e.
% what is left after we remove the impact from the
% covariates)
[~,res] = y_regress_ss(squeeze(V0(i,j,k,V0idx)),X);
[~,res_VX] = y_regress_ss(squeeze(V0_VX(i,j,k,V0idx)),X);
% The volumes following regression of covariates are created
V0Cov(i,j,k,V0idx)=res+mean(squeeze(V0(i,j,k,V0idx)));
V0Cov_VX(i,j,k,V0idx)=res_VX+mean(squeeze(V0_VX(i,j,k,V0idx)));
end
end
end
% Nan elements are replaced by 0 elements
WriteInformation_TB(fid,['Nan atlas: ',num2str(sum(sum(sum(sum(isnan(V0Cov)))))),'\t']);
WriteInformation_TB(fid,['Nanvoxel: ',num2str(sum(sum(sum(sum(isnan(V0Cov_VX)))))),'\n']);
V0Cov(isnan(V0Cov))=0;
V0Cov_VX(isnan(V0Cov_VX))=0;
%% 11. Writing NIFTI files for the regressed data (voxel-wise)
fprintf('Writing new volumes...');
% Will contain the data of interest (header)
mV0i=V0i_VX;
% For each volume over time, we creates the name of the header and
% write the output
for y=1:size(V0Cov_VX,4)
% Name creation
s=strfind(V0i_VX{y}.fname,'/');
last=s(size(s,2));
mV0i{y}.fname=strcat(V0i_VX{y}.fname(1:last),'Cov_',V0i_VX{y}.fname(last+1:end));
% Writing
spm_write_vol(mV0i{y}, V0Cov_VX(:,:,:,y));
fprintf('.');
end
%% 12. Atlas based averaging (smoothed data)
disp('Atlas-based averaging: AAL...');
[~, TCS_A1] = atlased_smoothing_TB(V0Cov(:,:,:,V0idx),Ai{1},A{1},...
codeBook{1},firstEV);
disp('Atlas-based averaging: Craddock 300...');
[~, TCS_A2] = atlased_smoothing_TB(V0Cov(:,:,:,V0idx),Ai{2},A{2},...
codeBook{2},firstEV);
disp('Atlas-based averaging: Craddock 900...');
[~, TCS_A3] = atlased_smoothing_TB(V0Cov(:,:,:,V0idx),Ai{3},A{3},...
codeBook{3},firstEV);
disp('Atlas-based averaging: James...');
[~, TCS_A4] = atlased_smoothing_TB(V0Cov(:,:,:,V0idx),Ai{4},A{4},...
codeBook{4},firstEV);
disp('Atlas-based averaging: Shen...');
[~, TCS_A5] = atlased_smoothing_TB(V0Cov(:,:,:,V0idx),Ai{5},A{5},...
codeBook{5},firstEV);
%% 13. Despiking
if Despike,
fprintf('Despiking\n');
for r=1:size(TCS1block,1)
tmp=TCS1block(r,:);
tmp(tmp > median(tmp) + 3*1.4826*mad(tmp,1)) = median(tmp) + 2.5*1.4826*mad(tmp,1);
tmp(tmp < median(tmp) - 3*1.4826*mad(tmp,1)) = median(tmp) - 2.5*1.4826*mad(tmp,1);
TCS1block(r,:)=tmp;
end
else
disp('Not despiking...');
end
%% 14. Bandpass filtering
if ~isempty(Band)
TCS1block = y_IdealFilter(TCS1block', TR, Band)'; % 1*2 array [low high], [0.01 0.1]
else
disp('Not filtering...');
end
disp('End of preprocessing !');
end
c4science · Help