reho.m

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

function [] = reho(NVolume,NVoxel, AMaskFilename)
% Calculate regional homogeneity (i.e. ReHo) from the 3D EPI images normalized by MNI template in SPM2. It is based on Matlab (5.3 or above) on windows 2000.
% FORMAT     function []   = reho(NVolume,NVoxel,bMask)
%                 NVolme   - The numbers of the volumes in the EPI images.
%                            (i.e. the length of the time series)
%                 NVoxel   - The number of the voxel for a given cluster
%                            during calculating the KCC (e.g. 27, 19, or 7);
%                            Recommand: NVoxel=27;
%                  bMask   - If the mask is used during calculating the KCC,
%                            bMmask = 1; if not, bMmask = 0. Here, four
%                            mask (reho_mask{02,03,12,13}.mat) were provided
%                            according to the normalized MNI image ('single_subj_T1.mnc')
%                            at SPM2. Recommand: bMask = 1. If you use other
%                            normalized items, please set bMask = 0.
% OUTPUT   Two analzye format file (rehomap.{hdr,img}) in the current directory.
% NOTE:    1.Please unzip reho.zip file and copy them to ...\matlab\spm2. 
%            After that, set ....\matlab\spm2 as a default path.
%          2 In the current directory, only the normalized EPI images in SPM2 that are
%            used to calculate the kcc are acceptable. Other "*.{hdr,img}" files should be 
%            removed from the directory!
%
% Example: reho(190,27,1);
% In the current directory, there are 190 volumes (i.e. NVolume = 190) EPI 
% functional images.During calculating the KCC, 27 voxel for a functinonal 
% cluster is chosen (NVoxel=27). Only the voxels within the brain are calculated,
% so bMask = 1. After that, you will get a analyze file (rehomap.{hdr,img}) in the
% current directory that denotes the regional homogeneity of the data set.

% Written by Yong He, April,2004
% Medical Imaging and Computing Group (MIC), National Laboratory of Pattern Recognition (NLPR),
% Chinese Academy of Sciences (CAS), China.
% E-mail: [e-mail]
% Copywrite (c) 2004, 

% Info on the approach based on the reho: 
% Zang YF, Jiang TZ, Lu YL, He Y and Tian LX, Regional Homogeneity Approach
% to fMRI Data Analysis. NeuroImage, Vol.22, No.1, 2004, 394-400.

% Info on the interesting and potential applications about the reho:
% He Y, Zang YF, Jiang TZ, Lu YL and Weng XC, Detection of Functional Networks
% in the Resting Brain. 2nd IEEE International Symposium on Biomedical Imaging:
% From Nano to Macro (ISBI'04), April 15-18, 2004, Arlington, USA. 

% Info on the above can be found in:
% http://www.nlpr.ia.ac.cn/english/mic/YongHe
% ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%Revised by Xiaowei Song, 20070421
%1.Add another parameter to allow user defined mask, so change some in code of selecting mask
%  user defined mask has priority 
% 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'
%2.Add waitbar for gui waiting and progress showing
%Revised by Xiaowei Song 20070903

if nargin < 1 error(' Input argument "NVolume,NVoxel,bMask" is undefined.'); end
if nargin > 3 error(' Error using ==> reho. Too many input arguments.'); end
theElapsedTime =cputime;
theOldWarnings =warning('off', 'all');
% judge the numbers of the volumes 
% -------------------------------------------------------------------------
tempfilename = dir;
count = 3; Number_volume = 0;
for count = 3:size(struct2cell(tempfilename),2)
    if tempfilename(count).name(end-3:end) == '.hdr' 
        if tempfilename(count).name(1:end-4) == tempfilename(count+1).name(1:end-4)
            Number_volume = Number_volume + 1;  
        else error('*.{hdr,img} should be pairwise. Please re-examin them.'); end
    end
end

% Examine the Nvoxel
% --------------------------------------------------------------------------
if NVoxel ~= 27 & NVoxel ~= 19 & NVoxel ~= 7 
    error('The second parameter should be 7, 19 or 27. Please re-exmamin it.');
end

%read the normalized functional images 
% -------------------------------------------------------------------------
fprintf('\n\t Read these 3D EPI functional images.\twait...');
[I,vsize,theImgFileList, Origin] =rest_to4d('.');
Long = size(I,4);

%mask selection, added by Xiaowei Song, 20070421
[pathstr, name, ext, versn] = fileparts(mfilename('fullpath'));
% examin the dimensions of the functional images and set mask 
M = size(I,1); N = size(I,2); O = size(I,3);
isize = [M N O];
mask=rest_loadmask(M, N, O, AMaskFilename);

%rank the 3d+time functional images voxel by voxel 
% -------------------------------------------------------------------------
fprintf('\n\t Calculate the rank of time series on voxel by voxel');
for ii = 1:M
    for jj = 1:N
        for kk = 1:O
			% 20071129 ~=0
            if mask(ii,jj,kk) ~=0 %== 255  % Long suggest adjustation 20070425,
				if Long<256,
					I(ii,jj,kk,:) = uint8(f_rank(squeeze(I(ii,jj,kk,:))));
				else
					I(ii,jj,kk,:) = uint16(f_rank(squeeze(I(ii,jj,kk,:))));
				end
            end							
        end %end k
    end %end j	
	rest_waitbar(ii/M, ...
			sprintf('Calculate the rank of time series\nwait...'), ...
			'ReHo Computing','Child','NeedCancelBtn');
	fprintf('.');
end %end i

% calulate the kcc for the data set
% ------------------------------------------------------------------------
fprintf('\t Calculate the kcc on voxel by voxel for the data set.\n');
K = zeros(M,N,O); 
switch NVoxel
    case 27  
        for i = 2:M-1
            for j = 2:N-1
                for k = 2:O-1                            
                    block = I(i-1:i+1,j-1:j+1,k-1:k+1,:);
                    mask_block = mask(i-1:i+1,j-1:j+1,k-1:k+1);
                    if min(min(min(mask_block))) > 0 
                        for ii = 1:Long
                            R_block(ii,:) = reshape(block(:,:,:,ii),1,27);
                        end
                        mask_R_block = R_block(:,reshape(mask_block,1,27) > 0);
                        K(i,j,k) = f_kendall(mask_R_block);  
                    end %end if	
                end%end k 
            end% end j			
			rest_waitbar(i/M, ...
					sprintf('Calculate the kcc\nwait...'), ...
					'ReHo Computing','Child','NeedCancelBtn');
        end%end i
        fprintf('\t The reho of the data set was finished.\n');
        fprintf('\t It has been written to the current directory (rehomap.{hdr, img}).\n');
        rest_writefile(K,'rehomap',isize,vsize,Origin, 'double');
    case 19  
        for i = 2:M-1
            for j = 2:N-1
                for k = 2:O-1                            
                    block = I(i-1:i+1,j-1:j+1,k-1:k+1,:);
                    mask_block = mask(i-1:i+1,j-1:j+1,k-1:k+1);
                    mask_block(1,1,1) = 2;    mask_block(1,3,1) = 2;    mask_block(3,1,1) = 2;    mask_block(3,3,1) = 2;  
                    mask_block(1,1,3) = 2;    mask_block(1,3,3) = 2;    mask_block(3,1,3) = 2;    mask_block(3,3,3) = 2;
                    if min(min(min(mask_block))) > 0  
                        for ii = 1:Long
                            R_block(ii,:) = reshape(block(:,:,:,ii),1,27);
                        end
                        mask_R_block = R_block(:,reshape(mask_block,1,27) > 2);
                        K(i,j,k) = f_kendall(mask_R_block);  
                    end%end if
                end%end k
            end%end j	
			rest_waitbar(i/M, ...
					sprintf('Calculate the kcc\nwait...'), ...
					'ReHo Computing','Child','NeedCancelBtn');
        end%end i
        fprintf('\t The reho of the data set was finished.\n');
        fprintf('\t It has been written to the current directory (rehomap.{hdr, img}).\n');
        rest_writefile(K,'rehomap',isize,vsize,Origin, 'double');
    case 7   
        for i = 2:M-1
            for j = 2:N-1
                for k = 2:O-1
                    block = I(i-1:i+1,j-1:j+1,k-1:k+1,:);
                    mask_block = mask(i-1:i+1,j-1:j+1,k-1:k+1);
                    mask_block(1,1,1) = 2;    mask_block(1,2,1) = 2;     mask_block(1,3,1) = 2;      mask_block(1,1,2) = 2;
                    mask_block(1,3,2) = 2;    mask_block(1,1,3) = 2;     mask_block(1,2,3) = 2;      mask_block(1,3,3) = 2; 
                    mask_block(2,1,1) = 2;    mask_block(2,3,1) = 2;     mask_block(2,1,3) = 2;      mask_block(2,3,3) = 2;   
                    mask_block(3,1,1) = 2;    mask_block(3,2,1) = 2;     mask_block(3,3,1) = 2;      mask_block(3,1,2) = 2;
                    mask_block(3,3,2) = 2;    mask_block(3,1,3) = 2;     mask_block(3,2,3) = 2;      mask_block(3,3,3) = 2;
                    if min(min(min(mask_block))) > 0 
                        for ii = 1:Long
                            R_block(ii,:) = reshape(block(:,:,:,ii),1,27);
                        end
                        mask_R_block = R_block(:,reshape(mask_block,1,27) > 2);
                        K(i,j,k) = f_kendall(mask_R_block);  
                    end%end if
                end%end k
            end%end j		
			rest_waitbar(i/M, ...
					sprintf('Calculate the kcc\nwait...'), ...
					'ReHo Computing','Child','NeedCancelBtn');
        end%end i
        fprintf('\t The reho of the data set was finished.\n');
        fprintf('\t It has been written to the current directory (rehomap.{hdr, img}).');
        rest_writefile(K,'rehomap',isize,vsize,Origin,'double');		
    otherwise
        error('The second parameter should be 7, 19 or 27. Please re-exmamin it.');
end %end switch
Ken = K;
theElapsedTime =cputime - theElapsedTime;
fprintf('\n\tRegional Homogeneity computation over, elapsed time: %g seconds\n', theElapsedTime);
warning(theOldWarnings);
end

% rank a time series
% -------------------------------------------------------------------------
function B = f_rank(A)
	SA = sort(A);
	for i = 1:length(A)
	    B(i,1) = mean(find(SA==A(i)));
	end
end
% calculate kcc for a time series
%---------------------------------------------------------------------------
function B = f_kendall(A)
	nk = size(A); n = nk(1); k = nk(2);
	SR = sum(A,2); SRBAR = mean(SR);
	S = sum(SR.^2) - n*SRBAR^2;
	B = 12*S/k^2/(n^3-n);
end

function Result = GetDirName(ADir)
	theDir =ADir;
	if strcmp(theDir(end),filesep)
		theDir=theDir(1:end-1);
	end	
	[tmp,Result]=fileparts(theDir);
end