fc

 Functional connectivity
AROIList would be treated as a mask in which time courses would be averaged to produce a new time course representing the ROI area
 Input:
     ADataDir            where the 3d+time dataset stay, and there should be 3d EPI functional image files. It must not contain / or \ at the end.
     AMaskFilename        the mask file name, I only compute the point within the mask
     AROIList        the mask list , ROI list definition
    AResultFilename        the output filename
    ACovariablesDef
 Output:
    AResultFilename    the filename of funtional connectivity result
-----------------------------------------------------------
    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
-----------------------------------------------------------
     Mail to Authors:  <a href="Dawnwei.Song@gmail.com">Xiaowei Song</a>; <a href="ycg.yan@gmail.com">Chaogan Yan</a> 
    Version=1.3;
    Release=20090321;
   Revised by YAN Chao-Gan, 080610. NIFTI compatible
   Last Revised by YAN Chao-Gan, 090321. Result data will be saved in the format 'single'.

0001 function [ResultMaps] = fc(ADataDir,AMaskFilename, AROIDef,AResultFilename, ACovariablesDef)
0002 % Functional connectivity
0003 %AROIList would be treated as a mask in which time courses would be averaged to produce a new time course representing the ROI area
0004 % Input:
0005 %     ADataDir            where the 3d+time dataset stay, and there should be 3d EPI functional image files. It must not contain / or \ at the end.
0006 %     AMaskFilename        the mask file name, I only compute the point within the mask
0007 %     AROIList        the mask list , ROI list definition
0008 %    AResultFilename        the output filename
0009 %    ACovariablesDef
0010 % Output:
0011 %    AResultFilename    the filename of funtional connectivity result
0012 %-----------------------------------------------------------
0013 %    Copyright(c) 2007~2010
0014 %    State Key Laboratory of Cognitive Neuroscience and Learning in Beijing Normal University
0015 %    Written by Xiao-Wei Song
0016 %    http://resting-fmri.sourceforge.net
0017 %-----------------------------------------------------------
0018 %     Mail to Authors:  <a href="Dawnwei.Song@gmail.com">Xiaowei Song</a>; <a href="ycg.yan@gmail.com">Chaogan Yan</a>
0019 %    Version=1.3;
0020 %    Release=20090321;
0021 %   Revised by YAN Chao-Gan, 080610. NIFTI compatible
0022 %   Last Revised by YAN Chao-Gan, 090321. Result data will be saved in the format 'single'.
0023 
0024     if nargin~=5,
0025         error(' Error using ==> fc. 4 arguments wanted.'); 
0026     end
0027         
0028     theElapsedTime =cputime;
0029     fprintf('\nComputing functional connectivity with:\t"%s"', ADataDir);
0030     [AllVolume,VoxelSize,theImgFileList, Header] =rest_to4d(ADataDir);
0031     % examin the dimensions of the functional images and set mask
0032     nDim1 = size(AllVolume,1); nDim2 = size(AllVolume,2); nDim3 = size(AllVolume,3);
0033     BrainSize = [nDim1 nDim2 nDim3]; 
0034     sampleLength = size(theImgFileList,1);
0035     
0036     %20070512    Saving a big 3D+time Dataset to small pieces by its first dimension to make this process run at least
0037     % put pieces of 4D dataset to the temp dir determined by the current time
0038     theTempDatasetDirName =sprintf('fc_%d', fix((1e4) *rem(now, 1) ));    
0039     theTempDatasetDir =[tempdir theTempDatasetDirName] ;
0040     ans=rmdir(theTempDatasetDir, 's');%suppress the error msg
0041     mkdir(tempdir, theTempDatasetDirName);    %Matlab 6.5 compatible
0042     
0043     AROIList =AROIDef;    
0044     if iscell(AROIDef),    %ROI wise, compute corelations between regions
0045         %ROI time course retrieving, 20070903
0046         theROITimeCourses =zeros(sampleLength, size(AROIDef,1));
0047         for x=1:size(AROIDef,1),
0048             fprintf('\n\t ROI time courses retrieving through "%s".', AROIDef{x});                
0049             IsDefinedROITimeCourse =0;
0050             if rest_SphereROI( 'IsBallDefinition', AROIDef{x})
0051                 %The ROI definition is a Ball definition
0052                 maskROI =rest_SphereROI( 'BallDefinition2Mask' , AROIDef{x}, BrainSize, VoxelSize, Header);
0053             elseif exist(AROIDef{x},'file')==2    % Make sure the Definition file exist
0054                 [pathstr, name, ext, versn] = fileparts(AROIDef{x});
0055                 if strcmpi(ext, '.txt'),
0056                     tmpX=load(AROIDef{x});
0057                     if size(tmpX,2)>1,
0058                         %Average all columns to make sure tmpX only contain one column
0059                         tmpX = mean(tmpX')';
0060                     end
0061                     theROITimeCourses(:, x) =tmpX;
0062                     IsDefinedROITimeCourse =1;
0063                 elseif strcmpi(ext, '.img')
0064                     %The ROI definition is a mask file
0065                     maskROI =rest_loadmask(nDim1, nDim2, nDim3, AROIDef{x});        
0066                 else
0067                     error(sprintf('REST doesn''t support the selected ROI definition now, Please check: \n%s', AROIDef{x}));
0068                 end
0069             else
0070                 error(sprintf('Wrong ROI definition, Please check: \n%s', AROIDef{x}));
0071             end
0072             
0073             if ~IsDefinedROITimeCourse,% I need retrieving the ROI averaged time course manualy
0074                 maskROI =find(maskROI);
0075                 % [rangeX, rangeY, rangeZ] = ind2sub(size(maskROI), find(maskROI));
0076                 % theTimeCourses =zeros(length(unique(rangeX)), length(unique(rangeY)), length(unique(rangeZ)));
0077                 for t=1:sampleLength,
0078                     theTimePoint = squeeze(AllVolume(:,:,:, t));
0079                     theTimePoint = theTimePoint(maskROI);
0080                     if ~isempty(theTimePoint),
0081                         theROITimeCourses(t, x) =mean(theTimePoint);
0082                     end
0083                 end    %The Averaged Time Course within the ROI now comes out! 20070903
0084             end%if ~IsDefinedROITimeCourse
0085         end%for
0086         %Save the ROI averaged time course to disk for further study
0087         [pathstr, name, ext, versn] = fileparts(AResultFilename);
0088         theROITimeCourseLogfile =[fullfile(pathstr,['ROI_', name]), '.txt'];
0089         save(theROITimeCourseLogfile, 'theROITimeCourses', '-ASCII', '-DOUBLE','-TABS')
0090         
0091         %If there are covariables
0092         theCovariables =[];
0093         if exist(ACovariablesDef.ort_file, 'file')==2,
0094             theCovariables =load(ACovariablesDef.ort_file);
0095             %Add polynomial in the baseline model according to 3dfim+.pdf
0096             thePolOrt=[];
0097             if ACovariablesDef.polort>=0,
0098                 thePolOrt =(1:sampleLength)';
0099                 thePolOrt =repmat(thePolOrt, [1, (1+ACovariablesDef.polort)]);
0100                 for x=1:(ACovariablesDef.polort+1),
0101                     thePolOrt(:, x) =thePolOrt(:, x).^(x-1) ;
0102                 end
0103             end
0104             
0105             %Regress out the covariables
0106             theCovariables =[thePolOrt, theCovariables];
0107             theROITimeCourses =Brain1D_RegressOutCovariables(theROITimeCourses, theCovariables);
0108             
0109         elseif ~isempty(ACovariablesDef.ort_file) && ~all(isspace(ACovariablesDef.ort_file)),
0110             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));
0111         end
0112         
0113         
0114         %Calcute the corelation matrix and save to a text file
0115         ResultMaps =corrcoef(theROITimeCourses);
0116         save([AResultFilename, '.txt'], 'ResultMaps', '-ASCII', '-DOUBLE','-TABS')
0117     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0118     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0119     else        %Voxel wise, compute corelations between one ROI time course and the whole brain
0120         %ROI time course retrieving, 20070903
0121         fprintf('\n\t ROI time course retrieving through "%s".', AROIDef);        
0122         AROITimeCourse = zeros(sampleLength, 1);
0123         IsDefinedROITimeCourse =0;
0124         if rest_SphereROI( 'IsBallDefinition', AROIDef)
0125             %The ROI definition is a Ball definition
0126             maskROI =rest_SphereROI( 'BallDefinition2Mask' , AROIDef, BrainSize, VoxelSize, Header);
0127         elseif exist(AROIDef,'file')==2    % Make sure the Definition file exist
0128             [pathstr, name, ext, versn] = fileparts(AROIDef);
0129             if strcmpi(ext, '.txt'),
0130                 tmpX=load(AROIDef);
0131                 if size(tmpX,2)>1,
0132                     %Average all columns to make sure tmpX only contain one column
0133                     tmpX = mean(tmpX')';
0134                 end
0135                 AROITimeCourse =tmpX;
0136                 IsDefinedROITimeCourse =1;
0137             elseif strcmpi(ext, '.img')
0138                 %The ROI definition is a mask file
0139                 maskROI =rest_loadmask(nDim1, nDim2, nDim3, AROIDef);        
0140             else
0141                 error(sprintf('REST doesn''t support the selected ROI definition now, Please check: \n%s', AROIDef));
0142             end
0143         else
0144             error(sprintf('Wrong ROI definition, Please check: \n%s', AROIDef));
0145         end
0146         if ~IsDefinedROITimeCourse,% I need retrieving the ROI averaged time course manualy
0147             maskROI =find(maskROI);
0148             % [rangeX, rangeY, rangeZ] = ind2sub(size(maskROI), find(maskROI));
0149             % theTimeCourses =zeros(length(unique(rangeX)), length(unique(rangeY)), length(unique(rangeZ)));
0150             for t=1:sampleLength,
0151                 theTimePoint = squeeze(AllVolume(:,:,:, t));
0152                 theTimePoint = theTimePoint(maskROI);
0153                 AROITimeCourse(t) =mean(theTimePoint);
0154             end    %The Averaged Time Course within the ROI now comes out! 20070903
0155             %Make sure the ROI averaged time course is an col vector
0156             AROITimeCourse =reshape(AROITimeCourse, sampleLength,1);
0157         end
0158         %Save the ROI averaged time course to disk for further study
0159         [pathstr, name, ext, versn] = fileparts(AResultFilename);
0160         theROITimeCourseLogfile =[fullfile(pathstr,['ROI_', name]), '.txt'];
0161         save(theROITimeCourseLogfile, 'AROITimeCourse', '-ASCII', '-DOUBLE','-TABS')
0162         
0163         %Save 3D+time Dataset's pieces to disk after ROI time course retrieved
0164         Save1stDimPieces(theTempDatasetDir, AllVolume, 'dim1_');
0165         clear AllVolume;%Free large memory
0166         
0167         %mask selection, added by Xiaowei Song, 20070421
0168         fprintf('\n\t Load mask "%s".', AMaskFilename);    
0169         mask=rest_loadmask(nDim1, nDim2, nDim3, AMaskFilename);
0170         
0171         fprintf('\n\t Build functional connectivity mask.\tWait...');
0172         mask =logical(mask);%Revise the mask to ensure that it contain only 0 and 1
0173         mask =    repmat(mask, [1, 1, 1, sampleLength]);    
0174         %Save mask pieces to disk to make this program at least run
0175         Save1stDimPieces(theTempDatasetDir, mask, 'mask_');    
0176         
0177         fprintf('\n\t Build Covariables time course.\tWait...');
0178         %If there are covariables, 20071002
0179         theCovariables =[];
0180         if exist(ACovariablesDef.ort_file, 'file')==2,
0181             theCovariables =load(ACovariablesDef.ort_file);
0182             %Add polynomial in the baseline model according to 3dfim+.pdf
0183             thePolOrt=[];
0184             if ACovariablesDef.polort>=0,
0185                 thePolOrt =(1:sampleLength)';
0186                 thePolOrt =repmat(thePolOrt, [1, (1+ACovariablesDef.polort)]);
0187                 for x=1:(ACovariablesDef.polort+1),
0188                     thePolOrt(:, x) =thePolOrt(:, x).^(x-1) ;
0189                 end
0190             end
0191             % switch ACovariablesDef.polort,
0192             % case 0,
0193                 % thePolOrt =ones(sampleLength, 1);
0194             % case 1,
0195                 % thePolOrt =[ones(sampleLength, 1), (1:sampleLength)'];
0196             % case 2,
0197                 % thePolOrt =[ones(sampleLength, 1), (1:sampleLength)', (1:sampleLength)'.^2];
0198             % otherwise
0199             % end
0200             %Regress out the covariables
0201             theCovariables =[thePolOrt, theCovariables];        
0202             
0203         elseif ~isempty(ACovariablesDef.ort_file) && ~all(isspace(ACovariablesDef.ort_file)),
0204             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));
0205         end
0206             
0207         fprintf('\n\t Functional connectivity is computing.\tWait...');        
0208         NumPieces_Dim1 =4;    %Constant number to divide the first dimension to "NumPieces_Dim1" pieces
0209         NumComputingCount =floor(nDim1/NumPieces_Dim1);
0210         if NumComputingCount< (nDim1/NumPieces_Dim1),
0211             NumComputingCount =NumComputingCount +1;
0212         else
0213         end
0214         for x=1:(NumComputingCount),    %20071129
0215         %for x=1:(floor(nDim1/NumPieces_Dim1) +1)
0216             rest_waitbar(x/(floor(nDim1/NumPieces_Dim1) +1), ...
0217                         'Computing Functional connectivity. Please wait...', ...
0218                         'REST working','Child','NeedCancelBtn');
0219                         
0220             %Load cached pieces of Datasets
0221             theFilename =fullfile(theTempDatasetDir, sprintf('dim1_%.8d', x));
0222             theDim1Volume4D =Load1stDimVolume(theFilename);
0223             theDim1Volume4D =double(theDim1Volume4D);
0224                     
0225             %Load and Apply the pieces' mask
0226             theFilename =fullfile(theTempDatasetDir, sprintf('mask_%.8d', x));
0227             theDim1Mask4D =Load1stDimVolume(theFilename);
0228             theDim1Volume4D(~theDim1Mask4D)=0;
0229 
0230             
0231             %I support multiple reference time course, and I will give each Ideals a Pearson correlation brain
0232             for y=1:size(AROITimeCourse, 2),
0233                 %Revise the dataset with covariables at first
0234                 if ~isempty(theCovariables),
0235                     AROITimeCourse(:, y) =Brain1D_RegressOutCovariables(AROITimeCourse(:, y), theCovariables);
0236                     theDim1Volume4D =Brain4D_RegressOutCovariables(theDim1Volume4D, theCovariables);
0237                 end
0238                 
0239                 %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
0240                 theCorrResult =PearsonCorrWithROI(AROITimeCourse(:, y), theDim1Volume4D);
0241                 
0242                 %Save to file
0243                 theFilename =fullfile(theTempDatasetDir, sprintf('result%.2d_%.8d', y, x));        
0244                 save(theFilename, 'theCorrResult');
0245             end
0246             fprintf('.');
0247         end
0248         clear theDim1Volume4D   theDim1Mask4D    theCorrResult;
0249         
0250         %Construct the 3D+time Dataset from files again
0251         fprintf('\n\t ReConstructing 3D Dataset Functional Connectivity.\tWait...');
0252         %Construct the correlation map's filenames, 20070905
0253         [pathstr, name, ext, versn] = fileparts(AResultFilename);
0254         ResultMaps =[];
0255         for  y=1:size(AROITimeCourse, 2),
0256             ResultMaps =[ResultMaps;{[pathstr, filesep ,name, sprintf('%.2d',y), ext]}];
0257         end
0258         %Reconstruct the Result correlation map from pieces
0259         for  y=1:size(AROITimeCourse, 2),
0260             theDataset3D=zeros(nDim1, nDim2, nDim3);
0261             for x=1:(NumComputingCount)
0262                 rest_waitbar(x/(floor(nDim1/NumPieces_Dim1)+1), ...
0263                             'Functional Connectivity 3D Brain reconstructing. Please wait...', ...
0264                             'REST working','Child','NeedCancelBtn');
0265                 
0266                 theFilename =fullfile(theTempDatasetDir, sprintf('result%.2d_%.8d', y, x));
0267                 %fprintf('\t%d',x);% Just for debugging
0268                 if x~=(floor(nDim1/NumPieces_Dim1)+1)
0269                     theDataset3D(((x-1)*NumPieces_Dim1+1):(x*NumPieces_Dim1),:,:)=Load1stDimVolume(theFilename);
0270                 else
0271                     theDataset3D(((x-1)*NumPieces_Dim1+1):end,:,:)=Load1stDimVolume(theFilename);
0272                 end        
0273                 fprintf('.');
0274             end
0275             
0276             if size(AROITimeCourse, 2)>1,
0277                 %Save every maps from result maps
0278                 fprintf('\n\t Saving Functional Connectivity maps: %s\tWait...', ResultMaps{y, 1});    
0279                 rest_writefile(single(theDataset3D), ...
0280                     ResultMaps{y, 1}, ...
0281                     BrainSize,VoxelSize,Header, 'single'); %Revised by YAN Chao-Gan, 090321. Result data will be stored in 'single' format. %'double');
0282             elseif size(AROITimeCourse, 2)==1,
0283                 %There will be no y loop, just one saving
0284                 %Save Functional Connectivity image to disk
0285                 fprintf('\n\t Saving Functional Connectivity map.\tWait...');    
0286                 rest_writefile(single(theDataset3D), ...
0287                     AResultFilename, ...
0288                     BrainSize,VoxelSize,Header, 'single'); %Revised by YAN Chao-Gan, 090321. Result data will be stored in 'single' format. %'double');
0289             end%end if
0290         end%end for
0291     end%voxel/region wise
0292     
0293 
0294     theElapsedTime =cputime - theElapsedTime;
0295     fprintf('\n\t Functional Connectivity compution over, elapsed time: %g seconds.\n', theElapsedTime);
0296     
0297 
0298     %After Band pass filter, remove the temporary files
0299     ans=rmdir(theTempDatasetDir, 's');%suppress the error msg
0300 %end
0301 
0302 %Save the 1st dimension of the 4D dataset to files
0303 function Save1stDimPieces(ATempDir, A4DVolume, AFilenamePrefix)
0304     NumPieces_Dim1 =4;    %Constant number to divide the first dimension to "NumPieces_Dim1" pieces
0305     NumComputingCount =floor(size(A4DVolume,1)/NumPieces_Dim1);
0306     if NumComputingCount< (size(A4DVolume,1)/NumPieces_Dim1),
0307         NumComputingCount =NumComputingCount +1;
0308     else
0309     end
0310     for x = 1:(NumComputingCount),
0311     %for x = 1:(floor(size(A4DVolume,1)/NumPieces_Dim1)+1)
0312         rest_waitbar(x/(floor(size(A4DVolume,1)/NumPieces_Dim1)+1), ...
0313                     'Cut one Big 3D+time Dataset into pieces of 3D+time Dataset Before Functional Connectivity. Please wait...', ...
0314                     'REST working','Child','NeedCancelBtn');
0315                     
0316         theFilename =fullfile(ATempDir, sprintf('%s%.8d',AFilenamePrefix, x));
0317         if x~=(floor(size(A4DVolume,1)/NumPieces_Dim1)+1)
0318             the1stDim = A4DVolume(((x-1)*NumPieces_Dim1+1):(x*NumPieces_Dim1), :,:,:);
0319         else
0320             the1stDim = A4DVolume(((x-1)*NumPieces_Dim1+1):end, :,:,:);
0321         end
0322         save(theFilename, 'the1stDim');         
0323     end    
0324 
0325 %Load the 1st dimension of the 4D dataset from files, return a Matrix not a struct
0326 function Result=Load1stDimVolume(AFilename)    
0327     Result =load(AFilename);
0328     theFieldnames=fieldnames(Result);    
0329     % Result =eval(sprintf('Result.%s',the1stField));%remove the struct variable to any named variable with a matrix
0330     Result = Result.(theFieldnames{1});
0331 
0332 %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
0333 %Dawnwei.Song, 20070903
0334 function ResultCorrCoef3DBrain=PearsonCorrWithROI(AROITimeCourse, ABrain4D);
0335     [nDim1, nDim2, nDim3, nDim4]=size(ABrain4D);
0336     %Remove the mean
0337     ABrain4D = ABrain4D - repmat(sum(ABrain4D,4)/nDim4,[1, 1, 1, nDim4]);
0338     AROITimeCourse =AROITimeCourse -repmat(sum(AROITimeCourse)/nDim4, [size(AROITimeCourse,1), size(AROITimeCourse,2)]);
0339     % For denominator
0340     ABrainSTD= squeeze(std(ABrain4D, 0, 4));
0341     ABrainSTD=ABrainSTD * std(AROITimeCourse);
0342     
0343     % (1*sampleLength) A matrix * B matrix (sampleLength * VoxelCount)
0344     ABrain4D =reshape(ABrain4D, nDim1*nDim2*nDim3, nDim4)';
0345     AROITimeCourse =reshape(AROITimeCourse, 1, nDim4);
0346     ABrain4D = AROITimeCourse * ABrain4D /(nDim4 -1);
0347     ABrain4D =squeeze(reshape(ABrain4D', nDim1, nDim2, nDim3, 1));
0348     
0349     % oldState =warning('query', 'MATLAB:divideByZero');
0350     % warning('off', 'MATLAB:divideByZero');
0351     ABrainSTD(find(ABrainSTD==0))=inf;%Suppress NaN to zero when denominator is zero
0352     ResultCorrCoef3DBrain =ABrain4D ./ABrainSTD;
0353     % ResultCorrCoef3DBrain(find(ABrainSTD==0)) =0; %Suppress NaN to zero when denominator is zero
0354     % warning(oldState.state, 'MATLAB:divideByZero');
0355     
0356 function Result =Brain4D_RegressOutCovariables(ABrain4D, ABasisFunc)
0357 %20070926, Regress some covariables out first
0358     %Result =( E - X(X'X)~X')Y
0359     [nDim1, nDim2, nDim3, nDim4]=size(ABrain4D);
0360     
0361     %Make sure the 1st dim of ABasisFunc is nDim4 long
0362     if size(ABasisFunc,1)~=nDim4, error('The length of Covariables don''t match with the volume.'); end
0363     
0364     % (1*sampleLength) A matrix * B matrix (sampleLength * VoxelCount)
0365     ABrain4D =reshape(ABrain4D, nDim1*nDim2*nDim3, nDim4)';
0366     Result =(eye(nDim4) - ABasisFunc * inv(ABasisFunc' * ABasisFunc)* ABasisFunc') * ABrain4D;
0367     %20071102 Bug fixed squeeze must not be excluded because nDim1 may be ONE !!!
0368     %Result =squeeze(reshape(Result', nDim1, nDim2, nDim3, nDim4));
0369     Result =reshape(Result', nDim1, nDim2, nDim3, nDim4);
0370 
0371 function Result =Brain1D_RegressOutCovariables(ABrain1D, ABasisFunc)
0372 %20070926, Regress some covariables out first
0373     %Result =( E - X(X'X)~X')Y
0374     %Make sure the input is a column vector
0375     % ABrain1D =reshape(ABrain1D, prod(size(ABrain1D)), 1);
0376     
0377     %Make sure the 1st dim of ABasisFunc is nDim4 long
0378     if size(ABasisFunc,1)~=length(ABrain1D), error('The length of Covariables don''t match with the volume.'); end
0379     
0380     % (1*sampleLength) A matrix * B matrix (sampleLength * VoxelCount)
0381     Result =(eye(size(ABrain1D, 1)) - ABasisFunc * inv(ABasisFunc' * ABasisFunc)* ABasisFunc') * ABrain1D;
0382