fc.m

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

			% case 1,
				% thePolOrt =[ones(sampleLength, 1), (1:sampleLength)'];
			% case 2,
				% thePolOrt =[ones(sampleLength, 1), (1:sampleLength)', (1:sampleLength)'.^2];
			% otherwise
			% end
			%Regress out the covariables
			theCovariables =[thePolOrt, theCovariables];		
			
		elseif ~isempty(ACovariablesDef.ort_file) && ~all(isspace(ACovariablesDef.ort_file)),
			warning(sprintf('\n\nCovariables definition text file "%s" doesn''t exist, please check! I won''t regress out the covariables this time.', ACovariablesDef.ort_file));
		end
			
		fprintf('\n\t Functional connectivity is computing.\tWait...');		
	    NumPieces_Dim1 =4;	%Constant number to divide the first dimension to "NumPieces_Dim1" pieces
		NumComputingCount =floor(nDim1/NumPieces_Dim1);
		if NumComputingCount< (nDim1/NumPieces_Dim1),
			NumComputingCount =NumComputingCount +1;
		else
		end
		for x=1:(NumComputingCount),	%20071129
		%for x=1:(floor(nDim1/NumPieces_Dim1) +1)
			rest_waitbar(x/(floor(nDim1/NumPieces_Dim1) +1), ...
						'Computing Functional connectivity. Please wait...', ...
						'REST working','Child','NeedCancelBtn');
						
			%Load cached pieces of Datasets
			theFilename =fullfile(theTempDatasetDir, sprintf('dim1_%.8d', x));
			theDim1Volume4D =Load1stDimVolume(theFilename);
			theDim1Volume4D =double(theDim1Volume4D);
					
			%Load and Apply the pieces' mask
			theFilename =fullfile(theTempDatasetDir, sprintf('mask_%.8d', x));
			theDim1Mask4D =Load1stDimVolume(theFilename);
			theDim1Volume4D(~theDim1Mask4D)=0;

			
			%I support multiple reference time course, and I will give each Ideals a Pearson correlation brain
			for y=1:size(AROITimeCourse, 2),
				%Revise the dataset with covariables at first
				if ~isempty(theCovariables),
					AROITimeCourse(:, y) =Brain1D_RegressOutCovariables(AROITimeCourse(:, y), theCovariables);
					theDim1Volume4D =Brain4D_RegressOutCovariables(theDim1Volume4D, theCovariables);
				end
				
				%Compute the Pearson coeffecient of ROI with a 3D+time Brain, return a 3D brain whose elements are the coefficients between the ROI averaged time course and the full 3D+time Brain
				theCorrResult =PearsonCorrWithROI(AROITimeCourse(:, y), theDim1Volume4D);
				
				%Save to file
				theFilename =fullfile(theTempDatasetDir, sprintf('result%.2d_%.8d', y, x));		
				save(theFilename, 'theCorrResult');
			end
			fprintf('.');
		end
		clear theDim1Volume4D   theDim1Mask4D	theCorrResult;
		
		%Construct the 3D+time Dataset from files again
		fprintf('\n\t ReConstructing 3D Dataset Functional Connectivity.\tWait...');
		%Construct the correlation map's filenames, 20070905
		[pathstr, name, ext, versn] = fileparts(AResultFilename);
		ResultMaps =[];
		for  y=1:size(AROITimeCourse, 2),
			ResultMaps =[ResultMaps;{[pathstr, filesep ,name, sprintf('%.2d',y), ext]}];
		end
		%Reconstruct the Result correlation map from pieces
		for  y=1:size(AROITimeCourse, 2),
			theDataset3D=zeros(nDim1, nDim2, nDim3);
			for x=1:(NumComputingCount)
				rest_waitbar(x/(floor(nDim1/NumPieces_Dim1)+1), ...
							'Functional Connectivity 3D Brain reconstructing. Please wait...', ...
							'REST working','Child','NeedCancelBtn');
				
				theFilename =fullfile(theTempDatasetDir, sprintf('result%.2d_%.8d', y, x));
				%fprintf('\t%d',x);% Just for debugging
				if x~=(floor(nDim1/NumPieces_Dim1)+1)
					theDataset3D(((x-1)*NumPieces_Dim1+1):(x*NumPieces_Dim1),:,:)=Load1stDimVolume(theFilename);
				else
					theDataset3D(((x-1)*NumPieces_Dim1+1):end,:,:)=Load1stDimVolume(theFilename);
				end		
				fprintf('.');
			end
			
			if size(AROITimeCourse, 2)>1,
				%Save every maps from result maps
				fprintf('\n\t Saving Functional Connectivity maps: %s\tWait...', ResultMaps{y, 1});	
				rest_writefile(theDataset3D, ...
					ResultMaps{y, 1}, ...
					BrainSize,VoxelSize,Origin, 'double');
			elseif size(AROITimeCourse, 2)==1,
				%There will be no y loop, just one saving
				%Save Functional Connectivity image to disk
				fprintf('\n\t Saving Functional Connectivity map.\tWait...');	
				rest_writefile(theDataset3D, ...
					AResultFilename, ...
					BrainSize,VoxelSize,Origin, 'double');
			end%end if
		end%end for	
	end%voxel/region wise
	

	theElapsedTime =cputime - theElapsedTime;
	fprintf('\n\t Functional Connectivity compution over, elapsed time: %g seconds.\n', theElapsedTime);
	

	%After Band pass filter, remove the temporary files
	ans=rmdir(theTempDatasetDir, 's');%suppress the error msg
%end

%Save the 1st dimension of the 4D dataset to files
function Save1stDimPieces(ATempDir, A4DVolume, AFilenamePrefix)
    NumPieces_Dim1 =4;	%Constant number to divide the first dimension to "NumPieces_Dim1" pieces
	NumComputingCount =floor(size(A4DVolume,1)/NumPieces_Dim1);
	if NumComputingCount< (size(A4DVolume,1)/NumPieces_Dim1),
		NumComputingCount =NumComputingCount +1;
	else
	end
	for x = 1:(NumComputingCount),
	%for x = 1:(floor(size(A4DVolume,1)/NumPieces_Dim1)+1)
		rest_waitbar(x/(floor(size(A4DVolume,1)/NumPieces_Dim1)+1), ...
					'Cut one Big 3D+time Dataset into pieces of 3D+time Dataset Before Functional Connectivity. Please wait...', ...
					'REST working','Child','NeedCancelBtn');
					
		theFilename =fullfile(ATempDir, sprintf('%s%.8d',AFilenamePrefix, x));
		if x~=(floor(size(A4DVolume,1)/NumPieces_Dim1)+1)
			the1stDim = A4DVolume(((x-1)*NumPieces_Dim1+1):(x*NumPieces_Dim1), :,:,:);
		else
			the1stDim = A4DVolume(((x-1)*NumPieces_Dim1+1):end, :,:,:);
		end
		save(theFilename, 'the1stDim'); 		
	end	

%Load the 1st dimension of the 4D dataset from files, return a Matrix not a struct
function Result=Load1stDimVolume(AFilename)	
	Result =load(AFilename);
	theFieldnames=fieldnames(Result);	
	% Result =eval(sprintf('Result.%s',the1stField));%remove the struct variable to any named variable with a matrix
	Result = Result.(theFieldnames{1});

%Compute the Pearson coeffecient of ROI with a 3D+time Brain, return a 3D brain whose elements are the coefficients between the ROI averaged time course and the full 3D+time Brain
%Dawnwei.Song, 20070903
function ResultCorrCoef3DBrain=PearsonCorrWithROI(AROITimeCourse, ABrain4D);
	[nDim1, nDim2, nDim3, nDim4]=size(ABrain4D);
	%Remove the mean
	ABrain4D = ABrain4D - repmat(sum(ABrain4D,4)/nDim4,[1, 1, 1, nDim4]);
	AROITimeCourse =AROITimeCourse -repmat(sum(AROITimeCourse)/nDim4, [size(AROITimeCourse,1), size(AROITimeCourse,2)]);
	% For denominator
	ABrainSTD= squeeze(std(ABrain4D, 0, 4));
	ABrainSTD=ABrainSTD * std(AROITimeCourse);
	
	% (1*sampleLength) A matrix * B matrix (sampleLength * VoxelCount)
	ABrain4D =reshape(ABrain4D, nDim1*nDim2*nDim3, nDim4)';
	AROITimeCourse =reshape(AROITimeCourse, 1, nDim4);
	ABrain4D = AROITimeCourse * ABrain4D /(nDim4 -1);
	ABrain4D =squeeze(reshape(ABrain4D', nDim1, nDim2, nDim3, 1));
	
	% oldState =warning('query', 'MATLAB:divideByZero');
	% warning('off', 'MATLAB:divideByZero');
	ABrainSTD(find(ABrainSTD==0))=inf;%Suppress NaN to zero when denominator is zero
	ResultCorrCoef3DBrain =ABrain4D ./ABrainSTD;
	% ResultCorrCoef3DBrain(find(ABrainSTD==0)) =0; %Suppress NaN to zero when denominator is zero
	% warning(oldState.state, 'MATLAB:divideByZero');
	
function Result =Brain4D_RegressOutCovariables(ABrain4D, ABasisFunc)
%20070926, Regress some covariables out first	
	%Result =( E - X(X'X)~X')Y
	[nDim1, nDim2, nDim3, nDim4]=size(ABrain4D);
	
	%Make sure the 1st dim of ABasisFunc is nDim4 long
	if size(ABasisFunc,1)~=nDim4, error('The length of Covariables don''t match with the volume.'); end
	
	% (1*sampleLength) A matrix * B matrix (sampleLength * VoxelCount)
	ABrain4D =reshape(ABrain4D, nDim1*nDim2*nDim3, nDim4)';
	Result =(eye(nDim4) - ABasisFunc * inv(ABasisFunc' * ABasisFunc)* ABasisFunc') * ABrain4D;
	%20071102 Bug fixed squeeze must not be excluded because nDim1 may be ONE !!!
	%Result =squeeze(reshape(Result', nDim1, nDim2, nDim3, nDim4));
	Result =reshape(Result', nDim1, nDim2, nDim3, nDim4);

function Result =Brain1D_RegressOutCovariables(ABrain1D, ABasisFunc)
%20070926, Regress some covariables out first	
	%Result =( E - X(X'X)~X')Y
	%Make sure the input is a column vector
	% ABrain1D =reshape(ABrain1D, prod(size(ABrain1D)), 1);
	
	%Make sure the 1st dim of ABasisFunc is nDim4 long
	if size(ABasisFunc,1)~=length(ABrain1D), error('The length of Covariables don''t match with the volume.'); end
	
	% (1*sampleLength) A matrix * B matrix (sampleLength * VoxelCount)	
	Result =(eye(size(ABrain1D, 1)) - ABasisFunc * inv(ABasisFunc' * ABasisFunc)* ABasisFunc') * ABrain1D;