Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Reduce volume memory usage #168

Merged
merged 18 commits into from
Mar 17, 2020
Merged

Reduce volume memory usage #168

Changes from 3 commits
Commits
Show all changes
18 commits
Select commit Hold shift + click to select a range
348c75d
Reduce volume memory usage
glasserm Mar 3, 2020
8ad2ca5
More logical and efficient sequence of steps for the the hp>0 condition
mharms Mar 4, 2020
b7303a8
Forgot to modify the fslmaths -bptf call to operate on the original f…
mharms Mar 4, 2020
e2a0f24
Change shell rm to matlab delete
coalsont Mar 4, 2020
9ac0eaa
change single(zeros(... to zeros(..., 'single')
coalsont Mar 4, 2020
97761a9
Add a couple linebreaks just to match other code
mharms Mar 5, 2020
624c8d1
MATLAB compile: functionhighpassandvariancenormalize
hodgem Mar 9, 2020
2430741
Update icaDim.m
glasserm Mar 10, 2020
bad46a5
Update icaDim.m
glasserm Mar 10, 2020
c05996a
Update icaDim.m
glasserm Mar 10, 2020
6edc725
Update icaDim.m
glasserm Mar 11, 2020
d122f4c
Some minor code restructuring to hopefully faciliate readability
mharms Mar 11, 2020
e0c8ca2
Decrease memory usage of functionhighpass by using 'range' instead of…
mharms Mar 12, 2020
a5b0cd2
Decrease memory usage of icaDim by using 'range' instead of 'std'; Ad…
mharms Mar 12, 2020
f1b2c5f
Delete 'data' once done with it
mharms Mar 12, 2020
9f08366
Fix some whitespace (tab vs space) misalignment
mharms Mar 12, 2020
ca736b1
Add reporting of percentage of masked voxels
mharms Mar 16, 2020
a75919f
MATLAB compile: functionhighpassandvariancenormalize
hodgem Mar 16, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
73 changes: 51 additions & 22 deletions ICAFIX/scripts/functionhighpassandvariancenormalize.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -93,13 +93,6 @@ function functionhighpassandvariancenormalize(TR,hp,fmri,WBC,varargin)
hpstring = ['_hp' pdstring num2str(hp)];
end

%% Load volume time series
if dovol > 0
cts=single(read_avw([fmri '.nii.gz']));
ctsX=size(cts,1); ctsY=size(cts,2); ctsZ=size(cts,3); ctsT=size(cts,4);
cts=reshape(cts,ctsX*ctsY*ctsZ,ctsT);
end

%% Load the motion confounds, and the CIFTI (if hp>=0) (don't need either if hp<0)
if hp>=0
confounds=load([fmri hpstring '.ica/mc/prefiltered_func_data_mcf.par']);
Expand All @@ -114,13 +107,28 @@ function functionhighpassandvariancenormalize(TR,hp,fmri,WBC,varargin)
%% Apply hp filtering of the motion confounds, volume (if requested), and CIFTI
if pdflag % polynomial detrend case
if dovol > 0
% Save and filter confounds, as a NIFTI, as expected by FIX
save_avw(reshape(confounds',size(confounds,2),1,1,size(confounds,1)),[fmri hpstring '.ica/mc/prefiltered_func_data_mcf_conf'],'f',[1 1 1 TR]);
confounds=detrendpoly(confounds,hp);
save_avw(reshape(confounds',size(confounds,2),1,1,size(confounds,1)),[fmri hpstring '.ica/mc/prefiltered_func_data_mcf_conf_hp'],'f',[1 1 1 TR]);


% Load volume time series and reduce to just the non-zero voxels (for memory efficiency)
ctsfull=single(read_avw([fmri '.nii.gz']));
ctsX=size(ctsfull,1); ctsY=size(ctsfull,2); ctsZ=size(ctsfull,3); ctsT=size(ctsfull,4);
ctsfull=reshape(ctsfull,ctsX*ctsY*ctsZ,ctsT);
ctsmask=std(ctsfull,[],2)>0;
cts=ctsfull(ctsmask,:);
clear ctsfull;

% Polynomial detrend
cts=detrendpoly(cts',hp)';
save_avw(reshape(cts,ctsX,ctsY,ctsZ,ctsT),[fmri hpstring '.nii.gz'],'f',[1 1 1 TR]);
% Use -d flag in fslcpgeom (even if not technically necessary) to reduce possibility of mistakes when editing script

% Write out result, restoring to original size
ctsfull=single(zeros(ctsX*ctsY*ctsZ,ctsT));
ctsfull(ctsmask,:)=cts;
save_avw(reshape(ctsfull,ctsX,ctsY,ctsZ,ctsT),[fmri hpstring '.nii.gz'],'f',[1 1 1 TR]);
clear ctsfull;
% Use -d flag in fslcpgeom (even if not technically necessary) to reduce possibility of mistakes when editing script
call_fsl(['fslcpgeom ' fmri '.nii.gz ' fmri hpstring '.nii.gz -d']);
end

Expand All @@ -129,26 +137,42 @@ function functionhighpassandvariancenormalize(TR,hp,fmri,WBC,varargin)

elseif hp>0 % "fslmaths -bptf" based filtering
if dovol > 0
% Save and filter confounds, as a NIFTI, as expected by FIX
save_avw(reshape(confounds',size(confounds,2),1,1,size(confounds,1)),[fmri hpstring '.ica/mc/prefiltered_func_data_mcf_conf'],'f',[1 1 1 TR]);
call_fsl(sprintf(['fslmaths ' fmri hpstring '.ica/mc/prefiltered_func_data_mcf_conf -bptf %f -1 ' fmri hpstring '.ica/mc/prefiltered_func_data_mcf_conf_hp'],0.5*hp/TR));

save_avw(reshape(cts,ctsX,ctsY,ctsZ,ctsT),[fmri hpstring '.nii.gz'],'f',[1 1 1 TR]);
call_fsl(['fslmaths ' fmri hpstring '.nii.gz -bptf ' num2str(0.5*hp/TR) ' -1 ' fmri hpstring '.nii.gz']);
cts=single(read_avw([fmri hpstring '.nii.gz']));
cts=reshape(cts,ctsX*ctsY*ctsZ,ctsT);
call_fsl(['fslcpgeom ' fmri '.nii.gz ' fmri hpstring '.nii.gz -d']);
% bptf filtering; no masking here, so this is probably memory inefficient
call_fsl(['fslmaths ' fmri '.nii.gz -bptf ' num2str(0.5*hp/TR) ' -1 ' fmri hpstring '.nii.gz']);
coalsont marked this conversation as resolved.
Show resolved Hide resolved
call_fsl(['fslcpgeom ' fmri '.nii.gz ' fmri hpstring '.nii.gz -d']);

% Load in result and reduce to just the non-zero voxels (for memory efficiency going forward)
ctsfull=single(read_avw([fmri hpstring '.nii.gz']));
ctsX=size(ctsfull,1); ctsY=size(ctsfull,2); ctsZ=size(ctsfull,3); ctsT=size(ctsfull,4);
ctsfull=reshape(ctsfull,ctsX*ctsY*ctsZ,ctsT);
ctsmask=std(ctsfull,[],2)>0;
cts=ctsfull(ctsmask,:);
clear ctsfull;
end

BOdimX=size(BO.cdata,1); BOdimZnew=ceil(BOdimX/100); BOdimT=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(single(read_avw('Atlas')),100*BOdimZnew,BOdimT); BO.cdata=grot(1:BOdimX,:); clear grot;
call_fsl('rm Atlas.nii.gz');
ciftisave(BO,[fmri '_Atlas' regstring hpstring '.dtseries.nii'],WBC);
call_fsl('rm Atlas.nii.gz');
coalsont marked this conversation as resolved.
Show resolved Hide resolved

elseif hp<0 % If no hp filtering, still need to at least demean the volumetric time series
if dovol > 0
cts=demean(cts')';

% Load volume time series and reduce to just the non-zero voxels (for memory efficiency)
ctsfull=single(read_avw([fmri '.nii.gz']));
ctsX=size(ctsfull,1); ctsY=size(ctsfull,2); ctsZ=size(ctsfull,3); ctsT=size(ctsfull,4);
ctsfull=reshape(ctsfull,ctsX*ctsY*ctsZ,ctsT);
ctsmask=std(ctsfull,[],2)>0;
cts=ctsfull(ctsmask,:);
clear ctsfull;

cts=demean(cts')';
end

end
Expand Down Expand Up @@ -176,14 +200,17 @@ function functionhighpassandvariancenormalize(TR,hp,fmri,WBC,varargin)
% (Only need these for the individual runs, which, per above, are expected to be passed in with hp requested).
if hp>=0
if dovol > 0
fname=[fmri hpstring '_vn.nii.gz'];
save_avw(reshape(Outcts.noise_unst_std,ctsX,ctsY,ctsZ,1),fname,'f',[1 1 1 1]);
call_fsl(['fslcpgeom ' fmri '_mean.nii.gz ' fname ' -d']);
fname=[fmri hpstring '_vn.nii.gz'];
vnfull=single(zeros(ctsX*ctsY*ctsZ,1));
vnfull(ctsmask)=Outcts.noise_unst_std;

save_avw(reshape(vnfull,ctsX,ctsY,ctsZ,1),fname,'f',[1 1 1 1]); clear vnfull;
call_fsl(['fslcpgeom ' fmri '_mean.nii.gz ' fname ' -d']);
end

VN=BO;
VN.cdata=OutBO.noise_unst_std;
fname=[fmri '_Atlas' regstring hpstring '_vn.dscalar.nii'];
fname=[fmri '_Atlas' regstring hpstring '_vn.dscalar.nii'];
disp(['saving ' fname]);
ciftisavereset(VN,fname,WBC);
end
Expand All @@ -195,7 +222,9 @@ function functionhighpassandvariancenormalize(TR,hp,fmri,WBC,varargin)
cts=cts./repmat(Outcts.noise_unst_std,1,ctsT);
% Use '_vnts' (volume normalized time series) as the suffix for the volumetric VN'ed TCS
fname=[fmri hpstring '_vnts.nii.gz'];
save_avw(reshape(cts,ctsX,ctsY,ctsZ,ctsT),fname,'f',[1 1 1 1]);
ctsfull=single(zeros(ctsX*ctsY*ctsZ,ctsT));
ctsfull(ctsmask,:)=cts;
save_avw(reshape(ctsfull,ctsX,ctsY,ctsZ,ctsT),fname,'f',[1 1 1 1]); clear ctsfull;
% N.B. Version of 'fslcpgeom' in FSL 6.0.0 requires a patch because it doesn't copy both the qform and sform faithfully
call_fsl(['fslcpgeom ' fmri '.nii.gz ' fname ' -d']);
end
Expand Down