reho

 Calculate regional homogeneity (i.e. ReHo) from the 3D EPI images.
 FORMAT     function []   = (ADataDir, NVoxel, AMaskFilename, AResultFilename)
 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.
   NVoxel              The number of the voxel for a given cluster during calculating the KCC (e.g. 27, 19, or 7); Recommand: NVoxel=27;
     AMaskFilename        the mask file name, I only compute the point within the mask
    AResultFilename        the output filename
 Output:
    AResultFilename    the filename of ReHo result

0001 function [] = reho(ADataDir, NVoxel, AMaskFilename, AResultFilename)
0002 % Calculate regional homogeneity (i.e. ReHo) from the 3D EPI images.
0003 % FORMAT     function []   = (ADataDir, NVoxel, AMaskFilename, AResultFilename)
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 %   NVoxel              The number of the voxel for a given cluster during calculating the KCC (e.g. 27, 19, or 7); Recommand: NVoxel=27;
0007 %     AMaskFilename        the mask file name, I only compute the point within the mask
0008 %    AResultFilename        the output filename
0009 % Output:
0010 %    AResultFilename    the filename of ReHo result
0011 
0012 % Written by Yong He, April,2004
0013 % Medical Imaging and Computing Group (MIC), National Laboratory of Pattern Recognition (NLPR),
0014 % Chinese Academy of Sciences (CAS), China.
0015 % E-mail: yhe@nlpr.ia.ac.cn
0016 % Copywrite (c) 2004,
0017 
0018 % Info on the approach based on the reho:
0019 % Zang YF, Jiang TZ, Lu YL, He Y and Tian LX, Regional Homogeneity Approach
0020 % to fMRI Data Analysis. NeuroImage, Vol.22, No.1, 2004, 394-400.
0021 
0022 % Info on the interesting and potential applications about the reho:
0023 % He Y, Zang YF, Jiang TZ, Lu YL and Weng XC, Detection of Functional Networks
0024 % in the Resting Brain. 2nd IEEE International Symposium on Biomedical Imaging:
0025 % From Nano to Macro (ISBI'04), April 15-18, 2004, Arlington, USA.
0026 
0027 % Info on the above can be found in:
0028 % http://www.nlpr.ia.ac.cn/english/mic/YongHe
0029 % ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
0030 %Revised by Xiaowei Song, 20070421
0031 %1.Add another parameter to allow user defined mask, so change some in code of selecting mask
0032 %  user defined mask has priority
0033 % And delete the old parameter 'bMask', if the AMaskFilename is null(bMask=0) or 'Default'(bMask=1), then I would use 'ones mask' or 'default mask'
0034 %2.Add waitbar for gui waiting and progress showing
0035 %Revised by Xiaowei Song 20070903
0036 % ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
0037 %   Revised by YAN Chao-Gan 080610: NIFTI compatible
0038 %   Revised by YAN Chao-Gan 080808: also support NIFTI images.
0039 %   Revised by YAN Chao-Gan, 090321. Result data will be saved in the format 'single'.
0040 %   Revised by YAN Chao-Gan, 090419. 1. Add Parameters: ADataDir and
0041 %        AResultFilename; Remove the old parameter 'NVolume'. 2. Fixed the bug
0042 %        of computing ReHo with 7 voxels or 19 voxels in a cluster.
0043 %   Algorithm re-written by YAN Chao-Gan, 090422. The algorithm has been re-written to speed up the calculation of ReHo.
0044 
0045 
0046 if nargin~=4
0047     error(' Error using ==> reho. 4 arguments wanted.');
0048 end
0049 theElapsedTime =cputime;
0050 
0051 % Examine the Nvoxel
0052 % --------------------------------------------------------------------------
0053 if NVoxel ~= 27 & NVoxel ~= 19 & NVoxel ~= 7 
0054     error('The second parameter should be 7, 19 or 27. Please re-exmamin it.');
0055 end
0056 
0057 %read the normalized functional images
0058 % -------------------------------------------------------------------------
0059 fprintf('\n\t Read these 3D EPI functional images.\twait...');
0060 [AllVolume,vsize,theImgFileList, Header] =rest_to4d(ADataDir);
0061 [nDim1 nDim2 nDim3 nDim4]=size(AllVolume);
0062 isize = [nDim1 nDim2 nDim3]; 
0063 mask=rest_loadmask(nDim1, nDim2, nDim3, AMaskFilename);
0064 
0065 %Algorithm re-written by YAN Chao-Gan, 090422. Speed up the calculation of ReHo.
0066 %Saving a big 3D+time Dataset to small pieces by its first dimension to make this process run at least
0067 % put pieces of 4D dataset to the temp dir determined by the current time
0068 theTempDatasetDirName =sprintf('ReHo_%d', fix((1e4) *rem(now, 1) ));
0069 theTempDatasetDir =[tempdir theTempDatasetDirName] ;
0070 ans=rmdir(theTempDatasetDir, 's');%suppress the error msg
0071 mkdir(tempdir, theTempDatasetDirName);    %Matlab 6.5 compatible
0072 
0073 Save1stDimPieces(theTempDatasetDir, AllVolume, 'dim1_');
0074 clear AllVolume;%Free large memory
0075 
0076 %rank the 3d+time functional images voxel by voxel
0077 % -------------------------------------------------------------------------
0078 fprintf('\n\t Calculate the rank of time series on voxel by voxel');
0079 
0080 NumPieces_Dim1=10;    %Constant number to divide the first dimension to "NumPieces_Dim1" pieces
0081 NumComputingCount =ceil(nDim1/NumPieces_Dim1);
0082 for x=1:(NumComputingCount),
0083     rest_waitbar(x/NumComputingCount, ...
0084         'Calculate the rank of time series on voxel by voxel. Please wait...', ...
0085         'REST working','Child','NeedCancelBtn');
0086 
0087     theFilename =fullfile(theTempDatasetDir, sprintf('dim1_%.8d', x));
0088     theDim1Volume4D =Load1stDimVolume(theFilename);
0089     theDim1Volume4D =double(theDim1Volume4D);
0090 
0091     % Algorithm re-written by YAN Chao-Gan, 090422. Speed up the calculation of ReHo.
0092     theDim1Volume4D=permute(theDim1Volume4D,[4,1,2,3]); % Change the Time Course to the first dimention
0093     [TimePoints,M,N,O]=size(theDim1Volume4D);
0094     [theDim1Volume4D,SortIndex] = sort(theDim1Volume4D,1);
0095     db=diff(theDim1Volume4D,1,1);
0096     clear theDim1Volume4D;
0097     db = db == 0;
0098     if x~=NumComputingCount
0099         mask_x = mask(((x-1)*NumPieces_Dim1+1):(x*NumPieces_Dim1), :,:);
0100         sumdb=squeeze(sum(db,1)).*mask_x;
0101     else
0102         mask_x = mask(((x-1)*NumPieces_Dim1+1):end, :,:);
0103         sumdb=sum(db,1);
0104         sumdb=reshape(sumdb,size(mask_x)).*mask_x;
0105     end
0106     SortedRanks=reshape(repmat([1:TimePoints]',M*N*O,1),[TimePoints,M,N,O]);
0107     if any(sumdb(:))
0108         TieAdjustIndex=find(sumdb);
0109         for i=1:length(TieAdjustIndex)
0110             [iM iN iO]=ind2sub([M N O],TieAdjustIndex(i));
0111             ranks=SortedRanks(:,iM,iN,iO);
0112             ties=db(:,iM,iN,iO);
0113             tieloc = [find(ties); TimePoints+2];
0114             maxTies = numel(tieloc);
0115             tiecount = 1;
0116             while (tiecount < maxTies)
0117                 tiestart = tieloc(tiecount);
0118                 ntied = 2;
0119                 while(tieloc(tiecount+1) == tieloc(tiecount)+1)
0120                     tiecount = tiecount+1;
0121                     ntied = ntied+1;
0122                 end
0123                 % Compute mean of tied ranks
0124                 ranks(tiestart:tiestart+ntied-1) = ...
0125                     sum(ranks(tiestart:tiestart+ntied-1)) / ntied;
0126                 tiecount = tiecount + 1;
0127             end
0128             SortedRanks(:,iM,iN,iO)=ranks;
0129         end
0130     end
0131     clear db sumdb;
0132     SortIndexBase=reshape(([1:M*N*O]'-1).*TimePoints,[M,N,O]);
0133     SortIndexBase=repmat(SortIndexBase,[1,1,1,TimePoints]);
0134     SortIndexBase=permute(SortIndexBase,[4,1,2,3]);
0135     SortIndex=SortIndex+SortIndexBase;
0136     clear SortIndexBase;
0137     theDim1Volume4D(SortIndex)=SortedRanks;
0138     theDim1Volume4D=reshape(theDim1Volume4D,size(SortedRanks));
0139     clear SortIndex SortedRanks;
0140     theDim1Volume4D=uint16(theDim1Volume4D); % Change to uint16 to get the same results of previous version.
0141     %Save to file
0142     theFilename =fullfile(theTempDatasetDir, sprintf('result_%.8d', x));
0143     save(theFilename, 'theDim1Volume4D');
0144     fprintf('.');
0145 
0146 end
0147 
0148 for x=1:(NumComputingCount)
0149     rest_waitbar(x/(floor(nDim1/NumPieces_Dim1)+1), ...
0150         'Rank 3D Brain reconstructing. Please wait...', ...
0151         'REST working','Child','NeedCancelBtn');
0152     theFilename =fullfile(theTempDatasetDir, sprintf('result_%.8d', x));
0153     if x==1
0154         I=Load1stDimVolume(theFilename);
0155     else
0156         I=cat(2,I,Load1stDimVolume(theFilename));
0157     end
0158     fprintf('.');
0159 end
0160 fprintf('\n');
0161 [TimePoints,M,N,O]=size(I);
0162 
0163 ans=rmdir(theTempDatasetDir, 's');%suppress the error msg
0164 
0165 % calulate the kcc for the data set
0166 % ------------------------------------------------------------------------
0167 fprintf('\t Calculate the kcc on voxel by voxel for the data set.\n');
0168 K = zeros(M,N,O); 
0169 switch NVoxel
0170     case 27  
0171         for i = 2:M-1
0172             for j = 2:N-1
0173                 for k = 2:O-1                            
0174                     block = I(:,i-1:i+1,j-1:j+1,k-1:k+1);
0175                     mask_block = mask(i-1:i+1,j-1:j+1,k-1:k+1);
0176                     if all(mask_block(:))
0177                         R_block=reshape(block,[],27); %Revised by YAN Chao-Gan, 090420. Speed up the calculation.
0178                         mask_R_block = R_block(:,reshape(mask_block,1,27) > 0);
0179                         K(i,j,k) = f_kendall(mask_R_block);  
0180                     end %end if
0181                 end%end k
0182             end% end j
0183             rest_waitbar(i/M, ...
0184                     sprintf('Calculate the kcc\nwait...'), ...
0185                     'ReHo Computing','Child','NeedCancelBtn');
0186         end%end i
0187         fprintf('\t The reho of the data set was finished.\n');
0188         rest_writefile(single(K),AResultFilename,isize,vsize,Header, 'single'); %Revised by YAN Chao-Gan, 090321. Result data will be stored in 'single' format. %'double');
0189     case 19  
0190         mask_cluster_19=ones(3,3,3);
0191         mask_cluster_19(1,1,1) = 0;    mask_cluster_19(1,3,1) = 0;    mask_cluster_19(3,1,1) = 0;    mask_cluster_19(3,3,1) = 0;
0192         mask_cluster_19(1,1,3) = 0;    mask_cluster_19(1,3,3) = 0;    mask_cluster_19(3,1,3) = 0;    mask_cluster_19(3,3,3) = 0;
0193         %Revised by YAN Chao-Gan, 090420. The element in the mask could be 1 other than 127. Fixed the bug of computing ReHo with 7 voxels or 19 voxels.
0194         for i = 2:M-1
0195             for j = 2:N-1
0196                 for k = 2:O-1                            
0197                     block = I(:,i-1:i+1,j-1:j+1,k-1:k+1);
0198                     mask_block = mask(i-1:i+1,j-1:j+1,k-1:k+1);
0199                     if all(mask_block(:))
0200                         mask_block=mask_block.*mask_cluster_19;
0201                         %Revised by YAN Chao-Gan, 090419. The element in the mask could be 1 other than 127. Fixed the bug of computing ReHo with 7 voxels or 19 voxels.
0202                         R_block=reshape(block,[],27);  %Revised by YAN Chao-Gan, 090420. Speed up the calculation.
0203                         mask_R_block = R_block(:,reshape(mask_block,1,27) > 0); %Revised by YAN Chao-Gan, 090419. The element in the mask could be 1 other than 127. Fixed the bug of computing ReHo with 7 voxels or 19 voxels. %> 2);
0204                         K(i,j,k) = f_kendall(mask_R_block);  
0205                     end%end if
0206                 end%end k
0207             end%end j
0208             rest_waitbar(i/M, ...
0209                     sprintf('Calculate the kcc\nwait...'), ...
0210                     'ReHo Computing','Child','NeedCancelBtn');
0211         end%end i
0212         fprintf('\t The reho of the data set was finished.\n');
0213         rest_writefile(single(K),AResultFilename,isize,vsize,Header, 'single'); %Revised by YAN Chao-Gan, 090321. Result data will be stored in 'single' format. %'double');
0214     case 7   
0215         mask_cluster_7=ones(3,3,3);
0216         mask_cluster_7(1,1,1) = 0;    mask_cluster_7(1,2,1) = 0;     mask_cluster_7(1,3,1) = 0;      mask_cluster_7(1,1,2) = 0;
0217         mask_cluster_7(1,3,2) = 0;    mask_cluster_7(1,1,3) = 0;     mask_cluster_7(1,2,3) = 0;      mask_cluster_7(1,3,3) = 0;
0218         mask_cluster_7(2,1,1) = 0;    mask_cluster_7(2,3,1) = 0;     mask_cluster_7(2,1,3) = 0;      mask_cluster_7(2,3,3) = 0;
0219         mask_cluster_7(3,1,1) = 0;    mask_cluster_7(3,2,1) = 0;     mask_cluster_7(3,3,1) = 0;      mask_cluster_7(3,1,2) = 0;
0220         mask_cluster_7(3,3,2) = 0;    mask_cluster_7(3,1,3) = 0;     mask_cluster_7(3,2,3) = 0;      mask_cluster_7(3,3,3) = 0;
0221         %Revised by YAN Chao-Gan, 090420. The element in the mask could be 1 other than 127. Fixed the bug of computing ReHo with 7 voxels or 19 voxels.
0222         for i = 2:M-1
0223             for j = 2:N-1
0224                 for k = 2:O-1
0225                     block = I(:,i-1:i+1,j-1:j+1,k-1:k+1);
0226                     mask_block = mask(i-1:i+1,j-1:j+1,k-1:k+1);
0227                     if all(mask_block(:))
0228                         mask_block=mask_block.*mask_cluster_7;
0229                         %Revised by YAN Chao-Gan, 090419. The element in the mask could be 1 other than 127. Fixed the bug of computing ReHo with 7 voxels or 19 voxels.
0230                         R_block=reshape(block,[],27); %Revised by YAN Chao-Gan, 090420. Speed up the calculation.
0231                         mask_R_block = R_block(:,reshape(mask_block,1,27) > 0); %Revised by YAN Chao-Gan, 090419. The element in the mask could be 1 other than 127. Fixed the bug of computing ReHo with 7 voxels or 19 voxels. %> 2);
0232                         K(i,j,k) = f_kendall(mask_R_block);  
0233                     end%end if
0234                 end%end k
0235             end%end j
0236             rest_waitbar(i/M, ...
0237                     sprintf('Calculate the kcc\nwait...'), ...
0238                     'ReHo Computing','Child','NeedCancelBtn');
0239         end%end i
0240         fprintf('\t The reho of the data set was finished.\n');
0241         rest_writefile(single(K),AResultFilename,isize,vsize,Header,'single'); %Revised by YAN Chao-Gan, 090321. Result data will be stored in 'single' format. %'double');
0242     otherwise
0243         error('The second parameter should be 7, 19 or 27. Please re-exmamin it.');
0244 end %end switch
0245 Ken = K;
0246 theElapsedTime =cputime - theElapsedTime;
0247 fprintf('\n\tRegional Homogeneity computation over, elapsed time: %g seconds\n', theElapsedTime);
0248 
0249 % calculate kcc for a time series
0250 %---------------------------------------------------------------------------
0251 function B = f_kendall(A)
0252 nk = size(A); n = nk(1); k = nk(2);
0253 SR = sum(A,2); SRBAR = mean(SR);
0254 S = sum(SR.^2) - n*SRBAR^2;
0255 B = 12*S/k^2/(n^3-n);
0256     
0257     
0258 function Save1stDimPieces(ATempDir, A4DVolume, AFilenamePrefix)
0259 %Save the 1st dimension of the 4D dataset to files
0260 NumPieces_Dim1=10;    %Constant number to divide the first dimension to "NumPieces_Dim1" pieces
0261 NumComputingCount =ceil(size(A4DVolume,1)/NumPieces_Dim1);
0262 for x = 1:(NumComputingCount),
0263     %for x = 1:(floor(size(A4DVolume,1)/NumPieces_Dim1)+1)
0264     rest_waitbar((x/NumComputingCount), ...
0265         'Cut one Big 3D+time Dataset into pieces of 3D+time Dataset Before ReHo. Please wait...', ...
0266         'REST working','Child','NeedCancelBtn');
0267 
0268     theFilename =fullfile(ATempDir, sprintf('%s%.8d',AFilenamePrefix, x));
0269     if x~=NumComputingCount
0270         the1stDim = A4DVolume(((x-1)*NumPieces_Dim1+1):(x*NumPieces_Dim1), :,:,:);
0271     else
0272         the1stDim = A4DVolume(((x-1)*NumPieces_Dim1+1):end, :,:,:);
0273     end
0274     save(theFilename, 'the1stDim');
0275 end
0276 
0277 function Result=Load1stDimVolume(AFilename)
0278 %Load the 1st dimension of the 4D dataset from files, return a Matrix not a struct
0279 Result =load(AFilename);
0280 theFieldnames=fieldnames(Result);
0281 % Result =eval(sprintf('Result.%s',the1stField));%remove the struct variable to any named variable with a matrix
0282 Result = Result.(theFieldnames{1});