rest_detrend.m

While resting-state fMRI is drawing more and more attention, there has not been a software for its d

function rest_detrend(ADataDir, APostfix)
%Removing the linear trend for REST by Xiao-Wei Song
%Usage: rest_detrend(ADataDir, APostfix)
%ADataDir			where the 3d+time dataset stay, and there should be 3d EPI functional image files. It must not contain / or \ at the end.
%------------------------------------------------------------------------------------------------------------------------------
% Remove linear trend 
% Save to ADataDir_APostfix	
%
%------------------------------------------------------------------------------------------------------------------------------
%	Copyright(c) 2007~2010
%	State Key Laboratory of Cognitive Neuroscience and Learning in Beijing Normal University
%	Written by Xiao-Wei Song 
%	http://resting-fmri.sourceforge.net
% [e-mail], Copyright 2007~2010
%------------------------------------------------------------------------------------------------------------------------------
% 	<a href="[e-mail]">Mail to Author</a>: Xiaowei Song
%	Version=1.0;
%	Release=20070903;

	fprintf('\nRemoving the linear trend:\n');
	[AllVolume,VoxelSize,theImgFileList, Origin] =rest_to4d(ADataDir);

	thePrecision ='int16';
    tmpData =double(squeeze(AllVolume(:, :, :, round(size(AllVolume, 4)/2))));
	if mean(abs(tmpData(0~=tmpData)))<100,	%I can't use mean because It use too much memory!
		thePrecision ='double';
	end
	
	% examin the dimensions of the functional images and set mask 
	nDim1 = size(AllVolume,1); nDim2 = size(AllVolume,2); nDim3 = size(AllVolume,3); nDim4 =size(AllVolume,4);
		
	% I have to create a for loop because detrend can support 2-dim at most
	for x=1:nDim1,
		rest_waitbar(x/nDim1, ...
				'Removing Linear Trend, wait...', ...
				'Detrend','Child','NeedCancelBtn');			
		
		oneAxialSlice =double(AllVolume(x, :, :, :));
		oneAxialSlice =reshape(oneAxialSlice, 1*nDim2*nDim3, nDim4)';
		oneAxialSlice =detrend(oneAxialSlice) +repmat(mean(oneAxialSlice), [size(oneAxialSlice,1), 1]);
		oneAxialSlice =reshape(oneAxialSlice', 1,nDim2,nDim3, nDim4);
		if strcmpi(thePrecision, 'int16'),
			AllVolume(x, :, :, :) =uint16(oneAxialSlice);
		else
			AllVolume(x, :, :, :) =(oneAxialSlice);
		end
		
		% dim3PlusTimeCourse =squeeze( AllVolume(x, :, :, :) );
		% theTimeCourse =double(dim3PlusTimeCourse'); %detrend only can do along the column, before detrend
		% dim3PlusTimeCourse =detrend(theTimeCourse);
		% dim3PlusTimeCourse =dim3PlusTimeCourse + ...
					% repmat(mean(theTimeCourse), [size(dim3PlusTimeCourse,1), 1]);
		% dim3PlusTimeCourse =dim3PlusTimeCourse'; %detrend only can do along the column, go back
		%% AllVolume(x, y, :, :) =uint16(dim3PlusTimeCourse);	%20071031, Dawnwei.Song revised!
		% if strcmpi(thePrecision, 'int16'),
			% AllVolume(x, :, :, :) =uint16(dim3PlusTimeCourse);
		% else
			% AllVolume(x, :, :, :) =(dim3PlusTimeCourse);
		% end
		
	end;

	if strcmp(ADataDir(end),filesep)==1,
		ADataDir=ADataDir(1:end-1);
	end	
	
	ADataDir =sprintf('%s%s',ADataDir,APostfix);
	ans=rmdir(ADataDir, 's');%suppress the error msg	
	[theParentDir,theOutputDirName]=fileparts(ADataDir);	
	mkdir(theParentDir, theOutputDirName);	%Matlab 6.5 compatible	
	
	sampleLength =size(theImgFileList,1);
	for x=1:sampleLength,
		rest_waitbar(x/sampleLength, ...
					sprintf('Saving to {hdr/img} pair files\nwait...'), ...
					'Remove linear trend Over','Child','NeedCancelBtn');
		rest_writefile(AllVolume(:, :, :, x), ...
			sprintf('%s%s%.8d', ADataDir, filesep,x), ...
			[nDim1,nDim2,nDim3],VoxelSize, Origin,'int16');
		if (mod(x,5)==0) %progress show
			fprintf('.');
		end
	end	  
	fprintf('\n');