rest_bandpass

Ideal Band pass filter for REST by Xiao-Wei Song
rest_bandpass(ADataDir, ...
 ASamplePeriod, ALowPass_HighCutoff, AHighPass_LowCutoff, ...
 ARetrend, ...
 ARemoveTrendBefore, ARemoveTrendAfter, ...                            
 AMaskFilename)
 Use Ideal rectangular filter to filter a 3d+time dataset
 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
     ALowPass_HighCutoff            low pass, high cutoff of the band, eg. 0.08
     AHighPass_LowCutoff            high pass,  low cutoff of the band, eg. 0.01
    ARetrend            'Yes' or 'No'.     if yes, then re add the linear trend after filtering after removing the trend
     AMaskFilename        the mask file name, compatible with old reho or reho_gui, can be 'Default' or 1, '' or 0, 'mask.mat', '../mask.img'
 Output:
     Create a new sibling-directory with ADataDir, and name as 'ADataDir_filtered', then put all filted images to the new sibling-directory
-----------------------------------------------------------
    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. Data in processing will not be converted to the format 'int16'.

0001 function [] = rest_bandpass(ADataDir, ...
0002                             ASamplePeriod, ALowPass_HighCutoff, AHighPass_LowCutoff, ...
0003                             ARetrend, ...
0004                             AMaskFilename)
0005 %Ideal Band pass filter for REST by Xiao-Wei Song
0006 %rest_bandpass(ADataDir, ...
0007                             % ASamplePeriod, ALowPass_HighCutoff, AHighPass_LowCutoff, ...
0008                             % ARetrend, ...
0009                             % ARemoveTrendBefore, ARemoveTrendAfter, ...
0010                             % AMaskFilename)
0011 % Use Ideal rectangular filter to filter a 3d+time dataset
0012 % Input:
0013 %     ADataDir            where the 3d+time dataset stay, and there should be 3d EPI functional image files. It must not contain / or \ at the end.
0014 %     ASamplePeriod        TR, or like the variable name
0015 %     ALowPass_HighCutoff            low pass, high cutoff of the band, eg. 0.08
0016 %     AHighPass_LowCutoff            high pass,  low cutoff of the band, eg. 0.01
0017 %    ARetrend            'Yes' or 'No'.     if yes, then re add the linear trend after filtering after removing the trend
0018 %     AMaskFilename        the mask file name, compatible with old reho or reho_gui, can be 'Default' or 1, '' or 0, 'mask.mat', '../mask.img'
0019 % Output:
0020 %     Create a new sibling-directory with ADataDir, and name as 'ADataDir_filtered', then put all filted images to the new sibling-directory
0021 %-----------------------------------------------------------
0022 %    Copyright(c) 2007~2010
0023 %    State Key Laboratory of Cognitive Neuroscience and Learning in Beijing Normal University
0024 %    Written by Xiao-Wei Song
0025 %    http://resting-fmri.sourceforge.net
0026 %-----------------------------------------------------------
0027 %     Mail to Authors:  <a href="Dawnwei.Song@gmail.com">Xiaowei Song</a>; <a href="ycg.yan@gmail.com">Chaogan Yan</a>
0028 %    Version=1.3;
0029 %    Release=20090321;
0030 %   Revised by YAN Chao-Gan 080610: NIFTI compatible
0031 %   Last Revised by YAN Chao-Gan, 090321. Data in processing will not be converted to the format 'int16'.
0032 
0033 
0034     if nargin~=6, help('rest_bandpass');error(' Error using ==> rest_bandpass. 6 arguments wanted.'); end
0035         
0036     tic;
0037     fprintf('\nIdeal rectangular filter:\t"%s"', ADataDir);
0038     [AllVolume,vsize,theImgFileList, Header] =rest_to4d(ADataDir);
0039 
0040     thePrecision ='double'; %Revised by YAN Chao-Gan, 090321. Data will not be converted to the format 'int16'. %thePrecision ='int16';
0041     tmpData =double(squeeze(AllVolume(:, :, :, round(size(AllVolume, 4)/2))));
0042     if mean(abs(tmpData(0~=tmpData)))<100,    %I can't use mean because It use too much memory!
0043         thePrecision ='double';
0044     end
0045     % examin the dimensions of the functional images and set mask
0046     nDim1 = size(AllVolume,1); nDim2 = size(AllVolume,2); nDim3 = size(AllVolume,3);
0047     isize = [nDim1 nDim2 nDim3]; 
0048         
0049     %20070512    Saving a big 3D+time Dataset to small pieces by its first dimension to make this process run at least
0050     % put pieces of 4D dataset to the temp dir determined by the current time
0051     theTempDatasetDirName =sprintf('BandPass_%d', fix((1e4) *rem(now, 1) ));    
0052     theTempDatasetDir =[tempdir theTempDatasetDirName] ;
0053     ans=rmdir(theTempDatasetDir, 's');%suppress the error msg
0054     mkdir(tempdir, theTempDatasetDirName);    %Matlab 6.5 compatible
0055     
0056     Save1stDimPieces(theTempDatasetDir, AllVolume, 'dim1_');
0057     clear AllVolume;%Free large memory
0058 
0059     %mask selection, added by Xiaowei Song, 20070421
0060     fprintf('\n\t Load mask "%s".', AMaskFilename);    
0061     mask=rest_loadmask(nDim1, nDim2, nDim3, AMaskFilename);
0062     
0063     fprintf('\n\t Build band pass filtered mask.\tWait...');
0064     sampleFreq      = 1/ASamplePeriod; 
0065     sampleLength = size(theImgFileList,1);
0066     paddedLength = rest_nextpow2_one35(sampleLength); %2^nextpow2(sampleLength);
0067     %paddedLength =sampleLength; %Don't pad any zero
0068     freqPrecision= sampleFreq/paddedLength;
0069     if rest_misc( 'GetMatlabVersion')<7.3
0070         warning off MATLAB:conversionToLogical;    
0071     end
0072     mask =logical(mask);%Revise the mask to ensure that it contain only 0 and 1
0073     if rest_misc( 'GetMatlabVersion')<7.3
0074         warning on MATLAB:conversionToLogical;
0075     end
0076     maskLowPass =    repmat(mask, [1, 1, 1, paddedLength]);
0077     maskHighPass=    maskLowPass;
0078     clear mask;
0079     %% GENERATE LOW PASS WINDOW    20070514, reference: fourior_filter.c in AFNI
0080     if (ALowPass_HighCutoff>=(sampleFreq/2))||(ALowPass_HighCutoff==0)        
0081         maskLowPass(:,:,:,:)=1;    %All pass
0082     elseif (ALowPass_HighCutoff>0)&&(ALowPass_HighCutoff< freqPrecision)        
0083         maskLowPass(:,:,:,:)=0;    % All stop
0084     else
0085         % Low pass, such as freq < 0.08 Hz
0086         idxCutoff    =round(ALowPass_HighCutoff *paddedLength *ASamplePeriod);
0087         idxCutoff2    =paddedLength+2 -idxCutoff;                %Center Index =(paddedLength/2 +1)
0088         %maskLowPass(:,:,:,1:idxCutoff)=1;            %Low pass, contain DC
0089         maskLowPass(:,:,:,idxCutoff+1:idxCutoff2-1)=0; %High eliminate
0090         %maskLowPass(:,:,:,idxCutoff2:paddedLength)=1;    %Low pass
0091     end
0092     %%GENERATE HIGH PASS WINDOW
0093     if (AHighPass_LowCutoff < freqPrecision)        
0094         maskHighPass(:,:,:,:)=1;    %All pass
0095     elseif (AHighPass_LowCutoff >= (sampleFreq/2))
0096         maskHighPass(:,:,:,:)=0;    % All stop
0097     else
0098         % high pass, such as freq > 0.01 Hz
0099         idxCutoff    =round(AHighPass_LowCutoff *paddedLength *ASamplePeriod);
0100         idxCutoff2    =paddedLength+2 -idxCutoff;                %Center Index =(paddedLength/2 +1)
0101         maskHighPass(:,:,:,1:idxCutoff-1)=0;    %Low eliminate
0102         %maskHighPass(:,:,:,idxCutoff:idxCutoff2)=1;    %High Pass
0103         maskHighPass(:,:,:,idxCutoff2+1:paddedLength)=0;    %Low eliminate
0104     end    
0105     %Combine the low pass mask and the high pass mask
0106     %maskLowPass(~maskHighPass)=0; %Don't combine because filter will not work when I only want low-pass or high-pass after combination, 20070517
0107     %Save mask pieces to disk to make this program at least run
0108     Save1stDimPieces(theTempDatasetDir, maskLowPass, 'fmLow_');    
0109     Save1stDimPieces(theTempDatasetDir, maskHighPass, 'fmHigh_');    
0110     clear maskLowPass maskHighPass; %Free large memory
0111     
0112     %20070513    remove trend --> FFT --> filter --> inverse FFT --> retrend
0113     if rest_misc( 'GetMatlabVersion')>=7.3
0114         fftw('dwisdom');        
0115     else        
0116     end    
0117     fprintf('\n\t Band Pass Filter working.\tWait...');        
0118     NumPieces_Dim1 =4;    %Constant number to divide the first dimension to "NumPieces_Dim1" pieces
0119     NumComputingCount =floor(nDim1/NumPieces_Dim1);
0120     if NumComputingCount< (nDim1/NumPieces_Dim1),
0121         NumComputingCount =NumComputingCount +1;
0122     else
0123     end
0124     for x=1:(NumComputingCount),
0125         rest_waitbar(x/(floor(nDim1/NumPieces_Dim1) +1), ...
0126                     'wait...', ...
0127                     'Filter is working','Child','NeedCancelBtn');
0128     
0129         %%Remove the linear trend first, ref: fourier_filter.c in AFNI, 20070509
0130         %Get every slope and intercept within the mask
0131         theFilename =fullfile(theTempDatasetDir, sprintf('dim1_%.8d', x));
0132         theDim1Volume4D =Load1stDimVolume(theFilename);
0133         theDim1Volume4D =double(theDim1Volume4D);
0134                 
0135         %Save the linear trend
0136         theTrend_Intercept=theDim1Volume4D(:,:,:, 1);
0137         theTrend_Slope= (theDim1Volume4D(:,:,:, end) -theTrend_Intercept) /double(sampleLength-1);
0138         for y=1:sampleLength
0139             %remove the linear trend first
0140             theDim1Volume4D(:,:,:, y)=theDim1Volume4D(:,:,:, y) -(theTrend_Intercept + y*theTrend_Slope);
0141         end
0142         theDim1Volume4D =cat(4,theDim1Volume4D,zeros(size(theDim1Volume4D,1),nDim2,nDim3,paddedLength -sampleLength));    %padded with zero
0143         
0144         %FFT
0145         theDim1Volume4D =fft(theDim1Volume4D, [], 4);
0146         %Low-pass Filter mask
0147         theFilename =fullfile(theTempDatasetDir, sprintf('fmLow_%.8d', x));
0148         theDim1FilterMask4D =Load1stDimVolume(theFilename);            
0149         %Apply the filter Low Pass
0150         theDim1Volume4D(~theDim1FilterMask4D)=0;
0151         
0152         %High-pass Filter mask
0153         theFilename =fullfile(theTempDatasetDir, sprintf('fmHigh_%.8d', x));
0154         theDim1FilterMask4D =Load1stDimVolume(theFilename);            
0155         %Apply the filter High Pass
0156         theDim1Volume4D(~theDim1FilterMask4D)=0;
0157         
0158         %inverse FFT
0159         theDim1Volume4D =ifft(theDim1Volume4D, [], 4);
0160         theDim1Volume4D =theDim1Volume4D(:,:,:, 1:sampleLength);%remove the padded parts
0161         %retrend the time course
0162         if strcmpi(ARetrend, 'Yes')
0163             for y=1:sampleLength
0164                 %add the linear trend after filter
0165                 theDim1Volume4D(:,:,:, y)=theDim1Volume4D(:,:,:, y) +double(theTrend_Intercept+y*theTrend_Slope);
0166             end
0167         end    
0168         %theDim1Volume4D =uint16(round(theDim1Volume4D));
0169         
0170         %Save to file
0171         theFilename =fullfile(theTempDatasetDir, sprintf('result_%.8d', x));        
0172         save(theFilename, 'theDim1Volume4D');         
0173         fprintf('.');
0174     end
0175     clear theDim1Volume4D theTrend_Intercept theTrend_Slope theDim1FilterMask4D;
0176     
0177     %Construct the 3D+time Dataset from files again
0178     fprintf('\n\t ReConstructing 3D+time Dataset.\tWait...');
0179     theDataset4D=zeros(nDim1, nDim2, nDim3, sampleLength);
0180     for x=1:(NumComputingCount)
0181         rest_waitbar(x/(floor(nDim1/NumPieces_Dim1)+1), ...
0182                     'wait...', ...
0183                     'After filter','Child','NeedCancelBtn');
0184                     
0185         theFilename =fullfile(theTempDatasetDir, sprintf('result_%.8d', x));
0186         %fprintf('\t%d',x);% Just for debugging
0187         if x~=(floor(nDim1/NumPieces_Dim1)+1)
0188             theDataset4D(((x-1)*NumPieces_Dim1+1):(x*NumPieces_Dim1),:,:,:)=Load1stDimVolume(theFilename);
0189         else
0190             theDataset4D(((x-1)*NumPieces_Dim1+1):end,:,:,:)=Load1stDimVolume(theFilename);
0191         end        
0192         fprintf('.');
0193     end
0194     
0195     %Save all images to disk
0196     fprintf('\n\t Saving filtered images.\tWait...');
0197     if strcmp(ADataDir(end),filesep)==1
0198         ADataDir=ADataDir(1:end-1);
0199     end    
0200     theResultOutputDir =sprintf('%s_filtered',ADataDir);
0201     ans=rmdir(theResultOutputDir, 's');%suppress the error msg
0202     [theParentDir,theOutputDirName]=fileparts(theResultOutputDir);    
0203     mkdir(theParentDir, theOutputDirName);    %Matlab 6.5 compatible
0204     
0205     %theDataset4D =uint16(theDataset4D);    %20071031 revised! Dawnwei.Song
0206     for x=1:sampleLength
0207         rest_waitbar(x/sampleLength, ...
0208                     sprintf('Saving to {hdr/img} pair files\nwait...'), ...
0209                     'Filter Over','Child','NeedCancelBtn');
0210         rest_writefile(single(theDataset4D(:, :, :, x)), ...
0211             sprintf('%s_filtered%s%.8d', ADataDir, filesep,x), ...
0212             isize,vsize,Header, 'single');  %Revised by YAN Chao-Gan, 090321. Filtered data will be stored in 'single' format. %thePrecision);
0213         if (mod(x,5)==0) %progress show
0214             fprintf('.');
0215         end
0216     end            
0217     fprintf('\n\t Band pass filter over.\n\t');
0218     toc;
0219 
0220     %After Band pass filter, remove the temporary files
0221     ans=rmdir(theTempDatasetDir, 's');%suppress the error msg
0222 
0223 %Save the 1st dimension of the 4D dataset to files
0224 function Save1stDimPieces(ATempDir, A4DVolume, AFilenamePrefix)
0225     NumPieces_Dim1 =4;    %Constant number to divide the first dimension to "NumPieces_Dim1" pieces
0226     NumComputingCount =floor(size(A4DVolume,1)/NumPieces_Dim1);
0227     if NumComputingCount< (size(A4DVolume,1)/NumPieces_Dim1),
0228         NumComputingCount =NumComputingCount +1;
0229     else
0230     end
0231     for x = 1:(NumComputingCount)
0232         rest_waitbar(x/(floor(size(A4DVolume,1)/NumPieces_Dim1)+1), ...
0233                     'wait...', ...
0234                     'Before Filter','Child','NeedCancelBtn');
0235         theFilename =fullfile(ATempDir, sprintf('%s%.8d',AFilenamePrefix, x));
0236         if x~=(floor(size(A4DVolume,1)/NumPieces_Dim1)+1)
0237             the1stDim = A4DVolume(((x-1)*NumPieces_Dim1+1):(x*NumPieces_Dim1), :,:,:);
0238         else
0239             the1stDim = A4DVolume(((x-1)*NumPieces_Dim1+1):end, :,:,:);
0240         end
0241         save(theFilename, 'the1stDim');         
0242     end    
0243     
0244 %Load the 1st dimension of the 4D dataset from files, return a Matrix not a struct
0245 function Result=Load1stDimVolume(AFilename)    
0246     Result =load(AFilename);
0247     theFieldnames=fieldnames(Result);    
0248     % Result =eval(sprintf('Result.%s',the1stField));%remove the struct variable to any named variable with a matrix
0249     Result = Result.(theFieldnames{1});
0250