alff

 Use ALFF method to compute the brain and return a ALFF brain map which reflects the "energy" of the voxels' BOLD signal
 FORMAT    function [] = alff(ADataDir,ASamplePeriod, ALowPass_HighCutoff, AHighPass_LowCutoff, 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.
     ASamplePeriod        TR, or like the variable name
     AHighPass_LowCutoff            the low edge of the pass band
     ALowPass_HighCutoff            the High edge of the pass band
     AMaskFilename        the mask file name, I only compute the point within the mask
    AResultFilename        the output filename
 Output:
    AResultFilename    the filename of ALFF 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 [] = alff(ADataDir,ASamplePeriod, ALowPass_HighCutoff, AHighPass_LowCutoff, AMaskFilename,AResultFilename)
0002 % Use ALFF method to compute the brain and return a ALFF brain map which reflects the "energy" of the voxels' BOLD signal
0003 % FORMAT    function [] = alff(ADataDir,ASamplePeriod, ALowPass_HighCutoff, AHighPass_LowCutoff, 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 %     ASamplePeriod        TR, or like the variable name
0007 %     AHighPass_LowCutoff            the low edge of the pass band
0008 %     ALowPass_HighCutoff            the High edge of the pass band
0009 %     AMaskFilename        the mask file name, I only compute the point within the mask
0010 %    AResultFilename        the output filename
0011 % Output:
0012 %    AResultFilename    the filename of ALFF result
0013 %-----------------------------------------------------------
0014 %    Copyright(c) 2007~2010
0015 %    State Key Laboratory of Cognitive Neuroscience and Learning in Beijing Normal University
0016 %    Written by Xiao-Wei Song
0017 %    http://resting-fmri.sourceforge.net
0018 %-----------------------------------------------------------
0019 %     Mail to Authors:  <a href="Dawnwei.Song@gmail.com">Xiaowei Song</a>; <a href="ycg.yan@gmail.com">Chaogan Yan</a>
0020 %    Version=1.3;
0021 %    Release=20090321;
0022 %   Revised by YAN Chao-Gan, 080610. NIFTI compatible
0023 %   Last Revised by YAN Chao-Gan, 090321. Result data will be saved in the format 'single'.
0024 
0025 
0026     if nargin~=6
0027         error(' Error using ==> alff. 6 arguments wanted.'); 
0028     end
0029         
0030     theElapsedTime =cputime;
0031     fprintf('\nComputing ALFF with:\t"%s"', ADataDir);
0032     [AllVolume,vsize,theImgFileList, Header] =rest_to4d(ADataDir);
0033     
0034     % examin the dimensions of the functional images and set mask
0035     [nDim1 nDim2 nDim3 nDim4]=size(AllVolume);
0036     %nDim1 = size(AllVolume,1); nDim2 = size(AllVolume,2); nDim3 = size(AllVolume,3);
0037     isize = [nDim1 nDim2 nDim3]; 
0038     
0039     %20070512    Saving a big 3D+time Dataset to small pieces by its first dimension to make this process run at least
0040     % put pieces of 4D dataset to the temp dir determined by the current time
0041     theTempDatasetDirName =sprintf('ALFF_%d', fix((1e4) *rem(now, 1) ));    
0042     theTempDatasetDir =[tempdir theTempDatasetDirName] ;
0043     ans=rmdir(theTempDatasetDir, 's');%suppress the error msg
0044     mkdir(tempdir, theTempDatasetDirName);    %Matlab 6.5 compatible
0045     
0046     Save1stDimPieces(theTempDatasetDir, AllVolume, 'dim1_');
0047     clear AllVolume;%Free large memory
0048 %mask selection, added by Xiaowei Song, 20070421
0049     fprintf('\n\t Load mask "%s".', AMaskFilename);    
0050     mask=rest_loadmask(nDim1, nDim2, nDim3, AMaskFilename);
0051     
0052     fprintf('\n\t Build ALFF filtered mask.\tWait...');
0053     sampleFreq      = 1/ASamplePeriod; 
0054     sampleLength = size(theImgFileList,1);
0055     paddedLength = rest_nextpow2_one35(sampleLength); %2^nextpow2(sampleLength);
0056     freqPrecision= sampleFreq/paddedLength;
0057         
0058     mask =logical(mask);%Revise the mask to ensure that it contain only 0 and 1
0059     maskLowPass =    repmat(mask, [1, 1, 1, paddedLength]);
0060     maskHighPass=    maskLowPass;
0061     clear mask;
0062     
0063     %20071226 ALFF parameters
0064     idxLowPass_HighCutoff    =round(ALowPass_HighCutoff *paddedLength *ASamplePeriod);
0065     %% GENERATE LOW PASS WINDOW    20070514, reference: fourior_filter.c in AFNI
0066     if (ALowPass_HighCutoff>=sampleFreq/2)||(ALowPass_HighCutoff==0)        
0067         maskLowPass(:,:,:,:)=1;    %All pass
0068     elseif (ALowPass_HighCutoff>0)&&(ALowPass_HighCutoff< freqPrecision)        
0069         maskLowPass(:,:,:,:)=0;    % All stop
0070     else
0071         % Low pass, such as freq < 0.08 Hz
0072         idxLowPass_HighCutoff    =round(ALowPass_HighCutoff *paddedLength *ASamplePeriod);
0073         idxLowPass_HighCutoff2    =paddedLength+2 -idxLowPass_HighCutoff;                %Center Index =(paddedLength/2 +1)
0074         %maskLowPass(:,:,:,1:idxCutoff)=1;            %Low pass, contain DC
0075         maskLowPass(:,:,:,idxLowPass_HighCutoff+1:idxLowPass_HighCutoff2-1)=0; %High eliminate
0076         %maskLowPass(:,:,:,idxCutoff2:paddedLength)=1;    %Low pass
0077     end
0078     
0079     %ALFF parameters 20071226
0080     idxHighPass_LowCutoff    =round(AHighPass_LowCutoff *paddedLength *ASamplePeriod);
0081     %%GENERATE HIGH PASS WINDOW
0082     if (AHighPass_LowCutoff < freqPrecision)    
0083         maskHighPass(:,:,:,:)=1;    %All pass
0084     elseif (AHighPass_LowCutoff >= sampleFreq/2)
0085         maskHighPass(:,:,:,:)=0;    % All stop
0086     else
0087         % high pass, such as freq > 0.01 Hz
0088         idxHighPass_LowCutoff    =round(AHighPass_LowCutoff *paddedLength *ASamplePeriod);
0089         idxHighPass_LowCutoff2    =paddedLength+2 -idxHighPass_LowCutoff;                %Center Index =(paddedLength/2 +1)
0090         maskHighPass(:,:,:,1:idxHighPass_LowCutoff-1)=0;    %Low eliminate
0091         %maskHighPass(:,:,:,idxCutoff:idxCutoff2)=1;    %High Pass
0092         maskHighPass(:,:,:,idxHighPass_LowCutoff2+1:paddedLength)=0;    %Low eliminate
0093     end    
0094     %Combine the low pass mask and the high pass mask
0095     %(~maskHighPass)=0;    %Don't combine because filter will not work when I only want low-pass or high-pass after combination, 20070517
0096     %Save mask pieces to disk to make this program at least run
0097     Save1stDimPieces(theTempDatasetDir, maskLowPass, 'fmLow_');    
0098     Save1stDimPieces(theTempDatasetDir, maskHighPass, 'fmHigh_');    
0099     clear maskLowPass maskHighPass; %Free large memory
0100 
0101     %20070513    remove trend --> FFT --> filter --> inverse FFT --> retrend
0102     if rest_misc( 'GetMatlabVersion')>=7.3
0103         fftw('dwisdom');
0104     end    
0105     fprintf('\n\t ALFF computing.\tWait...');        
0106     NumPieces_Dim1 =4;    %Constant number to divide the first dimension to "NumPieces_Dim1" pieces
0107     NumComputingCount =floor(nDim1/NumPieces_Dim1);
0108     if NumComputingCount< (nDim1/NumPieces_Dim1),
0109         NumComputingCount =NumComputingCount +1;
0110     else
0111     end
0112     for x=1:(NumComputingCount),
0113 %    for x=1:(floor(nDim1/NumPieces_Dim1) +1)
0114         rest_waitbar(x/(floor(nDim1/NumPieces_Dim1) +1), ...
0115                     'Computing ALFF. Please wait...', ...
0116                     'REST working','Child','NeedCancelBtn');
0117                     
0118         %%Remove the linear trend first, ref: fourier_filter.c in AFNI, 20070509
0119         %Get every slope and intercept within the mask
0120         theFilename =fullfile(theTempDatasetDir, sprintf('dim1_%.8d', x));
0121         theDim1Volume4D =Load1stDimVolume(theFilename);
0122         theDim1Volume4D =double(theDim1Volume4D);
0123                 
0124         %Save the linear trend
0125         % theTrend_Intercept=theDim1Volume4D(:,:,:, 1);
0126         % theTrend_Slope= (theDim1Volume4D(:,:,:, end) -theTrend_Intercept) /double(sampleLength-1);
0127         % for y=1:sampleLength % Does ALFF need this routine?
0128             % remove the linear trend first
0129             % theDim1Volume4D(:,:,:, y)=theDim1Volume4D(:,:,:, y) -(theTrend_Intercept + y*theTrend_Slope);
0130         % end
0131         
0132         %I must Detrend it first, 20070703
0133         % for xx=1:size(theDim1Volume4D,1), for yy=1:size(theDim1Volume4D,2),
0134             % dim3PlusTimeCourse =squeeze( theDim1Volume4D(xx, yy, :, :) );
0135             % detrend only support 2-D operations at most, so I have to do like this
0136             % dim3PlusTimeCourse =detrend(dim3PlusTimeCourse');
0137             % I didn't add the mean back this time, 20070703
0138             % theDim1Volume4D(xx, yy, :, :) =dim3PlusTimeCourse';
0139         % end;end;
0140         %20071110
0141         for xx=1:size(theDim1Volume4D,1),
0142             oneAxialSlice =double(theDim1Volume4D(xx, :, :, :));
0143             oneAxialSlice =reshape(oneAxialSlice, 1*nDim2*nDim3, nDim4)';
0144             oneAxialSlice =detrend(oneAxialSlice);% +repmat(mean(oneAxialSlice), [size(oneAxialSlice,1), 1]);
0145             oneAxialSlice =reshape(oneAxialSlice', 1,nDim2,nDim3, nDim4);
0146             theDim1Volume4D(xx, :, :, :) =(oneAxialSlice);
0147         end;
0148                 
0149         theDim1Volume4D =cat(4,theDim1Volume4D,zeros(size(theDim1Volume4D,1),nDim2,nDim3,paddedLength -sampleLength));    %padded with zero
0150         
0151         %FFT
0152         theDim1Volume4D =fft(theDim1Volume4D, [], 4);
0153         %Low-pass Filter mask
0154         theFilename =fullfile(theTempDatasetDir, sprintf('fmLow_%.8d', x));
0155         theDim1FilterMask4D =Load1stDimVolume(theFilename);            
0156         %Apply the filter Low Pass
0157         theDim1Volume4D(~theDim1FilterMask4D)=0;
0158         
0159         %High-pass Filter mask
0160         theFilename =fullfile(theTempDatasetDir, sprintf('fmHigh_%.8d', x));
0161         theDim1FilterMask4D =Load1stDimVolume(theFilename);            
0162         %Apply the filter High Pass
0163         theDim1Volume4D(~theDim1FilterMask4D)=0;
0164         
0165         %Get the amplitude only in one of the symmetric sides after FFT
0166         theDim1Volume4D =theDim1Volume4D(:, :, :, 2:(paddedLength/2+1));
0167         theDim1Volume4D =abs(theDim1Volume4D);
0168         %Get the Power Spectrum, double it because I only want one side of both sides
0169         theDim1Volume4D =2*(theDim1Volume4D .* theDim1Volume4D) /sampleLength;
0170         theDim1Volume4D(:,:,:, end) =theDim1Volume4D(:,:,:, end) /2;
0171         %The DC component didn't double because it didn't have its symetric side
0172         %theDim1Volume4D(:,:,:,1) =theDim1Volume4D(:,:,:,1)/2;
0173         %Get the Square root of the power spectrum between 0.01 and 0.08, i.e., ALFF
0174         theDim1Volume4D =sqrt(theDim1Volume4D);
0175         %Averaged ALFF across 0.01~0.08
0176         theDim1Volume4D =sum(theDim1Volume4D,4);
0177         theDim1Volume4D =theDim1Volume4D /(idxLowPass_HighCutoff -idxHighPass_LowCutoff +1);
0178         
0179         %Save to file
0180         theFilename =fullfile(theTempDatasetDir, sprintf('result_%.8d', x));        
0181         save(theFilename, 'theDim1Volume4D');         
0182         fprintf('.');
0183     end
0184     clear theDim1Volume4D theTrend_Intercept theTrend_Slope theDim1FilterMask4D oneAxialSlice;
0185     
0186     %Construct the 3D+time Dataset from files again
0187     fprintf('\n\t ReConstructing 3D Dataset ALFF.\tWait...');
0188     theDataset3D=zeros(nDim1, nDim2, nDim3);
0189     for x=1:(NumComputingCount)
0190         rest_waitbar(x/(floor(nDim1/NumPieces_Dim1)+1), ...
0191                     'ALFF 3D Brain reconstructing. Please wait...', ...
0192                     'REST working','Child','NeedCancelBtn');
0193         
0194         theFilename =fullfile(theTempDatasetDir, sprintf('result_%.8d', x));
0195         %fprintf('\t%d',x);% Just for debugging
0196         if x~=(floor(nDim1/NumPieces_Dim1)+1)
0197             theDataset3D(((x-1)*NumPieces_Dim1+1):(x*NumPieces_Dim1),:,:)=Load1stDimVolume(theFilename);
0198         else
0199             theDataset3D(((x-1)*NumPieces_Dim1+1):end,:,:)=Load1stDimVolume(theFilename);
0200         end        
0201         fprintf('.');
0202     end
0203     
0204     %Save ALFF image to disk
0205     fprintf('\n\t Saving ALFF map.\tWait...');    
0206     rest_writefile(single(theDataset3D), ...
0207         AResultFilename, ...
0208         isize,vsize,Header, 'single'); %Revised by YAN Chao-Gan, 090321. Result data will be stored in 'single' format. %'double');
0209 
0210     theElapsedTime =cputime - theElapsedTime;
0211     fprintf('\n\t ALFF compution over, elapsed time: %g seconds.\n', theElapsedTime);
0212 
0213     %After Band pass filter, remove the temporary files
0214     ans=rmdir(theTempDatasetDir, 's');%suppress the error msg
0215 %end
0216 
0217 %Save the 1st dimension of the 4D dataset to files
0218 function Save1stDimPieces(ATempDir, A4DVolume, AFilenamePrefix)
0219     NumPieces_Dim1 =4;    %Constant number to divide the first dimension to "NumPieces_Dim1" pieces
0220     NumComputingCount =floor(size(A4DVolume,1)/NumPieces_Dim1);
0221     if NumComputingCount< (size(A4DVolume,1)/NumPieces_Dim1),
0222         NumComputingCount =NumComputingCount +1;
0223     else
0224     end
0225     for x = 1:(NumComputingCount),
0226     %for x = 1:(floor(size(A4DVolume,1)/NumPieces_Dim1)+1)
0227         rest_waitbar(x/(floor(size(A4DVolume,1)/NumPieces_Dim1)+1), ...
0228                     'Cut one Big 3D+time Dataset into pieces of 3D+time Dataset Before ALFF. Please wait...', ...
0229                     'REST working','Child','NeedCancelBtn');
0230                     
0231         theFilename =fullfile(ATempDir, sprintf('%s%.8d',AFilenamePrefix, x));
0232         if x~=(floor(size(A4DVolume,1)/NumPieces_Dim1)+1)
0233             the1stDim = A4DVolume(((x-1)*NumPieces_Dim1+1):(x*NumPieces_Dim1), :,:,:);
0234         else
0235             the1stDim = A4DVolume(((x-1)*NumPieces_Dim1+1):end, :,:,:);
0236         end
0237         save(theFilename, 'the1stDim');         
0238     end    
0239 
0240 %Load the 1st dimension of the 4D dataset from files, return a Matrix not a struct
0241 function Result=Load1stDimVolume(AFilename)    
0242     Result =load(AFilename);
0243     theFieldnames=fieldnames(Result);    
0244     % Result =eval(sprintf('Result.%s',the1stField));%remove the struct variable to any named variable with a matrix
0245     Result = Result.(theFieldnames{1});