Click here to download this file

function [MeanWmDrift,Ci95WmDrift] = ExtractWhiteMatterDrift(WhiteTpm,fMask,fImgs)
% "Extract White Matter Drift"
% By Sam Berens (S.Berens@sussex.ac.uk).
%
% INPUTS:
%   - "WhiteTpm" : a strings denoting the file name of a white matter...
%   tissue probability map (TPM);
%   - "fMask" : a string denoting the file name of a binary mask for the...
%   functional data. This mask should be in registration with the white...
%   matter TPM;
%   - "fImgs" : a nx1 cell array containing strings denoting the file...
%   names of each scan in the functional time series;
%
% OUTPUTS:
%   - "MeanWmDrift" : an nx1 double coding for drift in the MR signal....
%   This is computed as the mean white matter intensity after voxel-wise...
%   z-scoring;
%   - "Ci95WmDrift" : an nx1 double indicating the 95% confidence...
%   interval of the MR signal drift;

% --- Load volumes:
Vol_WhiteTpm = spm_vol(WhiteTpm);
P_WhiteTpm = spm_imatrix(Vol_WhiteTpm.mat);
VoxDim_WhiteTpm = abs(P_WhiteTpm(7:9));

Vol_fMask = spm_vol(fMask);
P_fMask = spm_imatrix(Vol_fMask.mat);
VoxDim_fMask = abs(P_fMask(7:9));

RelativeSize = round(VoxDim_fMask ./ VoxDim_WhiteTpm);
xyzDiff = ((RelativeSize - 1) ./ 2) .* VoxDim_WhiteTpm;

% --- Produce a list of white matter coordinates for the functional images:
i_Coordinates_fVx = 0;
Coordinates_fVx = cell(i_Coordinates_fVx,1);
SamsDisp('Retrieving white matter coordinates in the functional volume');
for fVx_x = 1:Vol_fMask.dim(1,1)
    for fVx_y = 1:Vol_fMask.dim(1,2)
        for fVx_z = 1:Vol_fMask.dim(1,3)
            fSample = spm_sample_vol(Vol_fMask,fVx_x,fVx_y,fVx_z,0);
            if fSample > 0.9
                fCoordinate_MNI = (Vol_fMask.mat * [fVx_x,fVx_y,fVx_z,1]')';
                fCoordinate_MNI = fCoordinate_MNI(1,1:3);
                TpmMni_xStart = fCoordinate_MNI(1,1) - xyzDiff(1,1);
                TpmMni_xEnd = fCoordinate_MNI(1,1) + xyzDiff(1,1);
                TpmMni_yStart = fCoordinate_MNI(1,2) - xyzDiff(1,2);
                TpmMni_yEnd = fCoordinate_MNI(1,2) + xyzDiff(1,2);
                TpmMni_zStart = fCoordinate_MNI(1,3) - xyzDiff(1,3);
                TpmMni_zEnd = fCoordinate_MNI(1,3) + xyzDiff(1,3);
                i_TmpVoxelMatrix_x = 0;
                for TmpMni_x = TpmMni_xStart:VoxDim_WhiteTpm(1,1):TpmMni_xEnd
                    i_TmpVoxelMatrix_x = i_TmpVoxelMatrix_x + 1;
                    i_TmpVoxelMatrix_y = 0;
                    for TmpMni_y = TpmMni_yStart:VoxDim_WhiteTpm(1,2):TpmMni_yEnd
                        i_TmpVoxelMatrix_y = i_TmpVoxelMatrix_y + 1;
                        i_TmpVoxelMatrix_z = 0;
                        for TmpMni_z = TpmMni_zStart:VoxDim_WhiteTpm(1,3):TpmMni_zEnd
                            i_TmpVoxelMatrix_z = i_TmpVoxelMatrix_z + 1;
                            TmpCoordinate_Vx = (Vol_WhiteTpm.mat \ [TmpMni_x,TmpMni_y,TmpMni_z,1]')';
                            TmpCoordinate_Vx = round(TmpCoordinate_Vx(1,1:3));
                            Sample_Tpm = spm_sample_vol(Vol_WhiteTpm,TmpCoordinate_Vx(1,1),TmpCoordinate_Vx(1,2),TmpCoordinate_Vx(1,3),0);
                            TmpVoxelMatrix(i_TmpVoxelMatrix_x,i_TmpVoxelMatrix_y,i_TmpVoxelMatrix_z) =...
                                Sample_Tpm > 0.99; %#ok<AGROW> % Populate the TmpVoxelMatrix with bools.
                        end
                    end
                end
                if sum(sum(sum(TmpVoxelMatrix,1),2),3) == numel(TmpVoxelMatrix)
                    i_Coordinates_fVx = i_Coordinates_fVx + 1;
                    Coordinates_fVx{i_Coordinates_fVx,1} = [fVx_x,fVx_y,fVx_z];
                end
            end
        end
    end
end
n_Coordinates_fVx = i_Coordinates_fVx;
SamsDisp(sprintf('Identified %i unique white matter voxels in the functional volume',n_Coordinates_fVx));

% --- Populate "WmIntensityMatrix":
SamsDisp('Populating white matter intensity matrix:');
WmIntensityMatrix = zeros(size(fImgs,1),n_Coordinates_fVx);
dispstat('','init');
for i_fImgs = 1:size(fImgs,1)
    Vol_fImgs = spm_vol(fImgs{i_fImgs,1});
    for i_Coordinate = 1:size(Coordinates_fVx,1)
        WmIntensityMatrix(i_fImgs,i_Coordinate) =...
            spm_sample_vol(Vol_fImgs,Coordinates_fVx{i_Coordinate,1}(1,1),Coordinates_fVx{i_Coordinate,1}(1,2),Coordinates_fVx{i_Coordinate,1}(1,3),0);
    end
    dispstat(sprintf('                         ...Completed functional volume %i of %i',i_fImgs,size(fImgs,1)));
end

% --- Produce averages:
SamsDisp('Producing averages...');
MeanVoxelIntensities = nanmean(WmIntensityMatrix,1);
Wd = pwd;
cd(sprintf('%s%stoolbox%sstats%sstats',matlabroot,filesep,filesep,filesep));
SdVoxelIntensities = nanstd(WmIntensityMatrix,1,1);
WmDriftMatrix = zeros(size(WmIntensityMatrix));
for i_Scan = 1:1:size(WmIntensityMatrix,1)
    for i_Voxel = 1:1:size(WmIntensityMatrix,2)
        WmDriftMatrix(i_Scan,i_Voxel) = (WmIntensityMatrix(i_Scan,i_Voxel) - MeanVoxelIntensities(1,i_Voxel)) / SdVoxelIntensities(1,i_Voxel);
    end
end
MeanWmDrift = nanmean(WmDriftMatrix,2);
Ci95WmDrift = (nanstd(WmDriftMatrix,0,2)./sqrt(size(WmDriftMatrix,2))).*1.96;
cd(Wd);
SamsDisp('Done!');
return

function [] = SamsDisp( Str )

Time = clock;
Year = sprintf('%04d',Time(1,1));
Month = sprintf('%02d',Time(1,2));
Day = sprintf('%02d',Time(1,3));
Hour = sprintf('%02d',Time(1,4));
Min = sprintf('%02d',Time(1,5));
Sec = sprintf('%02d',round(Time(1,6)));

SrtToPrint = sprintf('[%s/%s/%s - %s:%s:%s]: %s',Day,Month,Year,Hour,Min,Sec,Str);
disp(SrtToPrint);

return

function dispstat(TXT,varargin)
% Prints overwritable message to the command line. If you dont want to keep
% this message, call dispstat function with option 'keepthis'. If you want to
% keep the previous message, use option 'keepprev'. First argument must be
% the message.
% IMPORTANT! In the firt call, option 'init' must be used for initialization purposes.
% Options:
%     'init'      this must be called in the begining. Otherwise, it can overwrite the previous outputs on the command line.
%     'keepthis'    the message will be persistent, wont be overwritable,
%     'keepprev'  the previous message wont be overwritten. New message will start from next line,
%     'timestamp' current time hh:mm:ss will be appended to the begining of the message.
% Example:
%   clc;
%   fprintf('12345677890\n');
%   dispstat('','init')      %Initialization. Does not print anything.
%   dispstat('Time stamp will be written over this text.'); % First output
%   dispstat('is current time.','timestamp','keepthis'); % Overwrites the previous output but this output wont be overwritten.
%   dispstat(sprintf('*********\nDeveloped by %s\n*********','Kasim')); % does not overwrites the previous output
%   dispstat('','timestamp','keepprev','keepthis'); % does not overwrites the previous output
%   dispstat('this wont be overwriten','keepthis');
%   dispstat('dummy dummy dummy');
%   dispstat('final stat');
% % Output:
%     12345677890
%     15:15:34 is current time.
%     *********
%     Developed by Kasim
%     *********
%     15:15:34
%     this wont be overwriten
%     final stat

% **********
% **** Options
keepthis = 0; % option for not overwriting
keepprev = 0;
timestamp = 0; % time stamp option
init = 0; % is it initialization step?
if ~isstr(TXT) %#ok<DISSTR>
    return
end
persistent prevCharCnt;
if isempty(prevCharCnt)
    prevCharCnt = 0;
end
if nargin == 0
    return
elseif nargin > 1
    for i = 2:nargin
        eval([varargin{i-1} '=1;']);
    end
end
if init == 1
    prevCharCnt = 0;
    return;
end
if isempty(TXT) && timestamp == 0
    return
end
if timestamp == 1
    c = clock; % [year month day hour minute seconds]
    txtTimeStamp = sprintf('%02d:%02d:%02d ',c(4),c(5),round(c(6)));
else
    txtTimeStamp = '';
end
if keepprev == 1
    prevCharCnt = 0;
end
% *************** Make safe for fprintf, replace control charachters
TXT = strrep(TXT,'%','%%');
TXT = strrep(TXT,'\','\\');
% *************** Print
TXT = [txtTimeStamp TXT '\n'];
fprintf([repmat('\b',1, prevCharCnt) TXT]);
nof_extra = length(strfind(TXT,'%%'));
nof_extra = nof_extra + length(strfind(TXT,'\\'));
nof_extra = nof_extra + length(strfind(TXT,'\n'));
prevCharCnt = length(TXT) - nof_extra; %-1 is for \n
if keepthis == 1
    prevCharCnt = 0;
end
return

Sam Berens