Home > rest_20090422 > f_alff.m

f_alff

PURPOSE ^

Use fALFF method to compute the brain and return a fALFF brain map.

SYNOPSIS ^

function [] = f_alff(ADataDir,ASamplePeriod, ALowPass_HighCutoff, AHighPass_LowCutoff, AMaskFilename,AResultFilename)

DESCRIPTION ^

 Use fALFF method to compute the brain and return a fALFF brain map.
 FORMAT    function [] = f_alff(ADataDir,ASamplePeriod, ALowPass_HighCutoff, AHighPass_LowCutoff, AMaskFilename,AResultFilename)
 Ref: Zou QH, Zhu CZ, Yang Y, Zuo XN, Long XY, Cao QJ, Wang YF, Zang YF (2008) An improved approach to detection of amplitude of low-frequency fluctuation (ALFF) for resting-state fMRI: fractional ALFF. Journal of neuroscience methods 172:137-141.
 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 fALFF 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. (1) Added the fALFF module, thank Dr. CHENG Wen-Lian for the helpful work. (2) Result data will be saved in the format 'single'.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

DOWNLOAD ^

f_alff.m

SOURCE CODE ^

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

Generated on Wed 29-Apr-2009 01:06:38 by m2html © 2005