Click here to download this file

function [] = CreateNativeEpiMask(In)
% "Create Native EPI Mask"
% By Sam Berens (S.Berens@sussex.ac.uk).
%
% INPUT: "In" should be a 1x1 structure with the following fields:
%   - "In.MeanEpi" : a strings denoting the file name of the mean EPI scan...
%   (or which ever EPI is being used as a reference) to base the mask on;
%   - "In.SaveName" : a string denoting the file name where the...
%   resultant mask should be saved;

if ~isfield(In,'MeanEpi')
    error('Reference image unspecified.');
end
if ~isfield(In,'SaveName')
    error('Output image filename unspecified.');
end

% --- 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;

% --- Get "VOL_MeanEpi":
Vol_MeanEpi = spm_vol(In.MeanEpi);

% --- Reshape data matrix and get histogram:
EpiVals = nan((Vol_MeanEpi.dim(1,1)*Vol_MeanEpi.dim(1,2)*Vol_MeanEpi.dim(1,3)),1);
for EpiVx_x = 1:Vol_MeanEpi.dim(1,1)
    for EpiVx_y = 1:Vol_MeanEpi.dim(1,2)
        for EpiVx_z = 1:Vol_MeanEpi.dim(1,3)
            Index = ((EpiVx_x-1)*Vol_MeanEpi.dim(1,2)*Vol_MeanEpi.dim(1,3))+((EpiVx_y-1)*Vol_MeanEpi.dim(1,3))+EpiVx_z;
            EpiVals(Index,1) = spm_sample_vol(Vol_MeanEpi,EpiVx_x,EpiVx_y,EpiVx_z,0);
        end
    end
end
[Density,Intensities] = ksdensity(EpiVals);

% --- 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 VOL and write to it:
Mat_NativeEpiMask = spm_read_vols(Vol_MeanEpi) > CutoffIntensity;
Vol_NativeEpiMask = Vol_MeanEpi;
if strcmp(In.SaveName(1,(end-3):end),'.img') || strcmp(In.SaveName(1,(end-3):end),'.nii')
    Vol_NativeEpiMask.fname = In.SaveName;
else
    Vol_NativeEpiMask.fname = sprintf('%s.img',In.SaveName);
end
Vol_NativeEpiMask.descrip = sprintf('sam - native epi mask (%f)',CutoffIntensity);
[~] = spm_write_vol(Vol_NativeEpiMask,Mat_NativeEpiMask);
return

Sam Berens