Click here to download this file

function [Vector_MeanAbsZs,Vector_AbsZsVar] = Sessa(In)
% "Sam's EPI Session Screening Algorithm (Sessa)"
% By Sam Berens (S.Berens@sussex.ac.uk); Thursday 31st July 2014.
%
% [Vector_MeanAbsZs,Vector_AbsZsVar] = Sessa(In)
%
% INPUT: "In" should be a 1x1 structure with the following fields:
%   - "In.Epis" : an nx1 cell array containing strings denoting the file...
%   names of n EPI scans to include in for screening;
%   - "In.SaveName" : [Optional] a string denoting the file name to which...
%   the resultant MATLAB figure should be saved;
%
% OUTPUTs: 2 nx1 vectors that can be used for diagnostic purposes (where...
% n = the number of EPIs):
%   - "Vector_MeanAbsZs" : Unsigned Z-scores averaged across non-empty EPI...
%   voxels. Values over 2 may indicate a problem;
%   - "Vector_AbsZsVar" : Unsigned standardised values of the variance...
%   across voxels within each volume. Values over 4 may indicate a problem;

if ~isfield(In,'Epis')
    error('Input volumes unspecified.');
end

% --- Return 4D matrix of EPIs (time along the 4th dimension):
FileNames_EPIs = In.Epis;
n_EPIs = size(FileNames_EPIs,1);
for i_EPI = 1:n_EPIs
    VOL_CurrentEpi = spm_vol(FileNames_EPIs{i_EPI,1});
    if i_EPI == 1
        Matrix_EpiData = nan(VOL_CurrentEpi.dim(1,1),VOL_CurrentEpi.dim(1,2),VOL_CurrentEpi.dim(1,3),n_EPIs);
    end
    Matrix_EpiData(:,:,:,i_EPI) = spm_read_vols(VOL_CurrentEpi);
    clear VOL_CurrentEpi;
end

% --- Produce mean EPI matrix, reshape and get histogram:
Matrix_MeanEpi = mean(Matrix_EpiData,4);
Vector_MeanEpi = reshape(Matrix_MeanEpi,(size(Matrix_MeanEpi,1)*size(Matrix_MeanEpi,2)*size(Matrix_MeanEpi,3)),1,1);
[Density,Intensities] = ksdensity(Vector_MeanEpi);

% --- Specify the filter specs:
Filter_Size = 20;
Filter_FWHM = 0.5;

% --- Produce the filter:
Filter_SD = (Filter_FWHM*Filter_Size) / (2*sqrt((2*log(2))));
Filter = zeros(Filter_Size,1);
Filter_Mean = Filter_Size/2;
for i = 1:size(Filter,1)
    Filter(i,1) = (1/(Filter_SD*(sqrt((2*pi))))) * exp(-1*(((i - Filter_Mean)^2)/(2*(Filter_SD^2))));
end
Filter_Max = max(Filter,[],1);
Filter = Filter ./ Filter_Max;

% --- Convolution:
SmoothedDensity = conv(Density,Filter,'same');

% --- Find Cut-off:
PassedPeak = 0;
CutoffIntensity = NaN;
for i_Bin = 1:(size(SmoothedDensity,2)-1)
    BinIndex = size(SmoothedDensity,2) - (i_Bin-1);
    if (~PassedPeak) && (SmoothedDensity(1,(BinIndex-1)) > SmoothedDensity(1,(BinIndex-0)))
        % Continue till find peak!
    elseif (~PassedPeak) && (SmoothedDensity(1,(BinIndex-1)) < SmoothedDensity(1,(BinIndex-0)))
        PassedPeak = 1;
    end
    if PassedPeak && (SmoothedDensity(1,(BinIndex-1)) > SmoothedDensity(1,(BinIndex-0)))
        CutoffIntensity = Intensities(1,BinIndex);
        break
    end
end

% --- Create new matrix and apply to 4D matrix of EPIs:
Matrix_EpiMask = Matrix_MeanEpi > CutoffIntensity;
for EpiVx_x = 1:size(Matrix_EpiMask,1)
    for EpiVx_y = 1:size(Matrix_EpiMask,2)
        for EpiVx_z = 1:size(Matrix_EpiMask,3)
            if ~Matrix_EpiMask(EpiVx_x,EpiVx_y,EpiVx_z)
                Matrix_EpiData(EpiVx_x,EpiVx_y,EpiVx_z,:) = nan(1,1,1,n_EPIs);
            end
        end
    end
end

% --- Calculate "Vector_MeanAbsZs":
Matrix_AbsZsEpiData = abs(zscore(Matrix_EpiData,1,4));
Vector_MeanAbsZs = nan(n_EPIs,1);
for i_EPI = 1:n_EPIs
    Vector_MeanAbsZs(i_EPI,1) = nanmean(squeeze(reshape(Matrix_AbsZsEpiData(:,:,:,i_EPI),(size(Matrix_AbsZsEpiData,1)*size(Matrix_AbsZsEpiData,2)*size(Matrix_AbsZsEpiData,3)),1,1,1)),1);
end

% --- Calculate "Vector_AbsZsVar":
Vector_Var = nan(n_EPIs,1);
for i_EPI = 1:n_EPIs
    Vector_Var(i_EPI,1) = nanvar(squeeze(reshape(Matrix_EpiData(:,:,:,i_EPI),(size(Matrix_EpiData,1)*size(Matrix_EpiData,2)*size(Matrix_EpiData,3)),1,1,1)),1,1);
end
Vector_AbsZsVar = abs(zscore(Vector_Var,1,1));

% --- Produce figures:
figure;
subplot(2,1,1);
plot(Vector_MeanAbsZs,'DisplayName','MeanAbsZs(LocalFx)','LineStyle','-','LineWidth',1,'Color','r');
legend({'MeanAbsZs(LocalFx)'},'Location','NorthEast');
xlabel('Volume Number');
ylabel('Statistic');
%
subplot(2,1,2);
plot(Vector_AbsZsVar,'DisplayName','AbsZsVar(GlobalFx)','LineStyle','-','LineWidth',1,'Color','b');
legend({'AbsZsVar(GlobalFx)'},'Location','NorthEast');
xlabel('Volume Number');
ylabel('Statistic');

% --- If requested, save plot:
if isfield(In,'SaveName')
    savefig(In.SaveName);
end
return

Sam Berens