function bmask = pm_brain_mask(P,flags)
% Calculate a brain mask
% FORMAT bmask = pm_brain_mask(P,flags)
%
% P - is a single pointer to a single image
%
% flags - structure containing various options
%         fwhm     - fwhm of smoothing kernel for generating mask
%         nerode   - number of erosions
%         thresh   - threshold for smoothed mask
%         ndilate  - number of dilations
%
%__________________________________________________________________________
%
% Inputs
% A single *.img conforming to SPM data format (see 'Data Format').
%
% Outputs
% Brain mask in a matrix
%__________________________________________________________________________
%
% The brain mask is generated by segmenting the image into GM, WM and CSF,
% adding these components together then thresholding above zero.
% A morphological opening is performed to get rid of stuff left outside of
% the brain. Any leftover holes are filled.
%__________________________________________________________________________

% Chloe Hutton
% Copyright (C) 2008-2022 Wellcome Centre for Human Neuroimaging


if nargin < 2 || isempty(flags)
%  flags.template=fullfile(spm('Dir'),'toolbox','FieldMap','T1.nii');
   flags.fwhm     = 5;
   flags.nerode   = 2;
   flags.ndilate  = 4;
   flags.thresh   = 0.5;
%  flags.reg      = 0.02;
%  flags.graphics = 0;
end

disp('Segmenting and extracting brain...');
%seg_flags.estimate.reg=flags.reg;
%seg_flags.graphics = flags.graphics;

% Updated to use renamed version of spm_segment
%VO    = pm_segment(P.fname,flags.template,seg_flags);
%bmask = double(VO(1).dat)+double(VO(2).dat)+double(VO(3).dat)>0;
bmask = get_bmask(P.fname);
bmask = open_it(bmask,flags.nerode,flags.ndilate); % Do opening to get rid of scalp

% Calculate kernel in voxels:
vxs   = sqrt(sum(P.mat(1:3,1:3).^2));
fwhm  = repmat(flags.fwhm,1,3)./vxs;
bmask = fill_it(bmask,fwhm,flags.thresh); % Do fill to fill holes

OP         = P;
OP.fname   = spm_file(P.fname,'prefix','bmask');
OP.descrip = sprintf('Mask:erode=%d,dilate=%d,fwhm=%d,thresh=%1.1f',flags.nerode,flags.ndilate,flags.fwhm,flags.thresh);
spm_write_vol(OP,bmask);
%__________________________________________________________________________

%__________________________________________________________________________
function bm = get_bmask(fnam)
% Use newer segmentation (ie less than 20 years old) to segment
tpm = fullfile(spm('dir'),'tpm','TPM.nii');
job = struct(...
    'channel', struct('vols',{{fnam}},'biasreg',0.001,'biasfwhm',60),...
    'tissue',  struct('tpm',cellfun(@(n)[tpm ',' num2str(n)],...
                                    num2cell(1:6),'UniformOutput',false),...
                      'ngaus',{1,1,2,3,4,2},...
                      'native',{[1 0],[1 0],[1 0],[0 0],[0 0],[0 0]},...
                      'warped',[0 0]),...
    'warp',    struct('mrf',1,'cleanup',1,'reg',[0 0.001 0.5 0.05 0.2],...
                      'affreg','mni','fwhm',0,'samp',4,'vox',NaN,'bb',nan(2,3)));
c   = spm_preproc_run(job,'inmem'); % Note that output is uint8
bm  = (c{1}{1} + c{1}{2} + c{1}{3}) > 255*0.25;
%__________________________________________________________________________

%__________________________________________________________________________
function vol = open_it(vol,ne,nd)
% Do a morphological opening. This consists of an erosion, followed by
% finding the largest connected component, followed by a dilation.

% Do an erosion then a connected components then a dilation
% to get rid of stuff outside brain.
for i=1:ne
   vol = spm_erode(double(vol));
end
vol = connect_it(vol);
for i=1:nd
   vol = spm_dilate(double(vol));
end
%__________________________________________________________________________

%__________________________________________________________________________
function ovol = fill_it(vol,k,thresh)
% Do morpholigical fill. This consists of finding the largest connected
% component and assuming that is outside of the head. All the other
% components are set to 1 (in the mask). The result is then smoothed by k
% and thresholded by thresh.
ovol = vol;

% Need to find connected components of negative volume
vol = ~vol;
[vol,NUM] = spm_bwlabel(double(vol),26);

% Now get biggest component and assume this is outside head..
pnc   = 0;
maxnc = 1;
for i=1:NUM
   nc = size(find(vol==i),1);
   if nc>pnc
      maxnc = i;
      pnc   = nc;
   end
end

% We know maxnc is largest cluster outside brain, so lets make all the
% others = 1.
for i=1:NUM
    if i~=maxnc
       ovol(vol==i) = 1;
    end
end

spm_smooth(ovol,ovol,k);
ovol = ovol>thresh;

%_______________________________________________________________________

%_______________________________________________________________________
function ovol = connect_it(vol)
% Find connected components and return the largest one.

[vol,NUM] = spm_bwlabel(double(vol),26);

% Get biggest component
pnc   = 0;
maxnc = 1;
for i=1:NUM
   nc=size(find(vol==i),1);
   if nc>pnc
      maxnc = i;
      pnc   = nc;
   end
end
ovol = (vol==maxnc);
