hmri_create_pm_brain_mask.m

function bmask = hmri_create_pm_brain_mask(P,flags)
% Calculate a brain mask in hMRI Toolbox
% This is pm_brain_mask (SPM12/toolbox/FieldMap) modified for the hMRI
% toolbox. Calls hmri_create_pm_segment instead of pm_segment. This is the
% only modification, syntax is unchanged otherwise.
%
% FORMAT bmask = hmri_create_pm_brain_mask(P,flags)
%
% P - is a single pointer to a single image
%
% flags - structure containing various options
%         template - which template for segmentation
%         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. 
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Chloe Hutton
% $Id: pm_brain_mask.m 5014 2012-10-24 10:56:28Z guillaume $

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...');

% % In pm_brain_mask, was:
% seg_flags.estimate.reg=flags.reg;
% seg_flags.graphics = flags.graphics;
% VO=pm_segment(P.fname,flags.template,seg_flags);
% % In hmri_create_pm_brain_mask, replaced by:
VO = hmri_create_pm_segment(P.fname);

bmask=double(VO(1).dat)+double(VO(2).dat)+double(VO(3).dat)>0;

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 ovol=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
   nvol=spm_erode(double(vol));
   vol=nvol;
end
nvol=connect_it(vol);
vol=nvol;
for i=1:nd
   nvol=spm_dilate(double(vol));
   vol=nvol;
end

ovol=nvol;
%__________________________________________________________________________

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