Added MH changes for fix_3_clean RAM footprint

CHANGELOG.txt

+17/03/20	1.06.15
+*		Lowered RAM footprint for fix_3_clean
 14/02/20	1.06.14
 *		Generate denoised motion regressors

fix_3_clean.m

@@ -68,17 +68,26 @@ if exist('Atlas.dtseries.nii','file') == 2
     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]);
+    BO.cdata = [];
     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; BO.cdata=BO.cdata+repmat(meanBO,1,size(BO.cdata,2));
+    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
+%%%%  read NIFTI version of the data, reducing to just the non-zero voxels (for memory efficiency)
 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)';
+  ctsfull=read_avw('filtered_func_data');
+  ctsX=size(ctsfull,1); ctsY=size(ctsfull,2); ctsZ=size(ctsfull,3); ctsT=size(ctsfull,4); 
+  % Note: reshape is a "memory-free" operation, but transpose isn't. So wait to transpose until after masking.
+  ctsfull=reshape(ctsfull,ctsX*ctsY*ctsZ,ctsT);
+  % Note: Use 'range' to identify non-zero voxels (which is very memory efficient)
+  % rather than 'std' (which requires additional memory equal to the size of the input)
+  ctsmask = range(ctsfull, 2) > 0;
+  cts=ctsfull(ctsmask,:)';  % Note: after transpose, cts has dimensions of [time space]
+  clear ctsfull;
 end
 
 %%%%  read and prepare motion confounds
@@ -95,6 +104,7 @@ if aggressive == 1
   sprintf('aggressive cleanup')
   confounds=[confounds ICA(:,DDremove)];
   if DOvol
+    % This regression peaks at a memory of > 2x (seemingly ~ 2.5x) the size of cts
     cts = cts - (confounds * (pinv(confounds,1e-6) * cts));
   end
   if DObrainord
@@ -124,7 +134,10 @@ end
 
 %%%% save cleaned data to file
 if DOvol
-  save_avw(reshape(cts',ctsX,ctsY,ctsZ,ctsT),'filtered_func_data_clean','f',[1 1 1 1]);
+  ctsfull=zeros(ctsX*ctsY*ctsZ,ctsT,class(cts));
+  ctsfull(ctsmask,:)=cts';
+  save_avw(reshape(ctsfull,ctsX,ctsY,ctsZ,ctsT),'filtered_func_data_clean','f',[1 1 1 1]);
+  clear ctsfull;
   call_fsl('fslcpgeom filtered_func_data filtered_func_data_clean');
 end
 if DObrainord
@@ -137,7 +150,9 @@ 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]);
+  vnfull=zeros(ctsX*ctsY*ctsZ,1,class(cts));
+  vnfull(ctsmask)=cts';
+  save_avw(reshape(vnfull,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

fixversion

-1.06.14
+1.06.15