function [segmentation] = ft_datatype_segmentation(segmentation, varargin)

% FT_DATATYPE_SEGMENTATION describes the FieldTrip MATLAB structure for segmented

% voxel-based data and atlasses. A segmentation can either be indexed or

% probabilistic (see below).

%

% A segmentation is a volumetric description which is usually derived from an

% anatomical MRI, which describes for each voxel the tissue type. It for example

% distinguishes between white matter, grey matter, csf, skull and skin. It is mainly

% used for masking in visualization, construction of volume conduction models and for

% construction of cortical sheets. An volume-based atlas is basically a very detailed

% segmentation with an anatomical label for each voxel.

%

% For example, the AFNI TTatlas+tlrc segmented brain atlas (which can be created

% with FT_READ_ATLAS) looks like this

%

% dim: [161 191 141] the size of the 3D volume in voxels

% transform: [4x4 double] affine transformation matrix for mapping the voxel coordinates to head coordinate system

% coordsys: 'tal' the transformation matrix maps the voxels into this (head) coordinate system

% unit: 'mm' the units in which the coordinate system is expressed

% brick0: [161x191x141 uint8] integer values from 1 to N, the value 0 means unknown

% brick1: [161x191x141 uint8] integer values from 1 to M, the value 0 means unknown

% brick0label: {Nx1 cell}

% brick1label: {Mx1 cell}

%

% An example segmentation with binary values that can be used for construction of a

% BEM volume conduction model of the head looks like this

%

% dim: [256 256 256] the dimensionality of the 3D volume

% transform: [4x4 double] affine transformation matrix for mapping the voxel coordinates to head coordinate system

% coordsys: 'ctf' the transformation matrix maps the voxels into this (head) coordinate system

% unit: 'mm' the units in which the coordinate system is expressed

% brain: [256x256x256 logical] binary map representing the voxels which belong to the brain

% skull: [256x256x256 logical] binary map representing the voxels which belong to the skull

% scalp: [256x256x256 logical] binary map representing the voxels which belong to the scalp

%

% An example of a whole-brain anatomical MRI that was segmented using FT_VOLUMESEGMENT

% looks like this

%

% dim: [256 256 256] the size of the 3D volume in voxels

% transform: [4x4 double] affine transformation matrix for mapping the voxel coordinates to head coordinate system

% coordsys: 'ctf' the transformation matrix maps the voxels into this (head) coordinate system

% unit: 'mm' the units in which the coordinate system is expressed

% gray: [256x256x256 double] probabilistic map of the gray matter

% white: [256x256x256 double] probabilistic map of the white matter

% csf: [256x256x256 double] probabilistic map of the cerebrospinal fluid

%

% The examples above demonstrate that a segmentation can be indexed, i.e. consisting

% of subsequent integer numbers (1, 2, ...) or probabilistic, consisting of real

% numbers ranging from 0 to 1 that represent probabilities between 0% and 100%. An

% extreme case is one where the probability is either 0 or 1, in which case the

% probability can be represented as a binary or logical array.

%

% The only difference to the volume data representation is that the segmentation

% structure contains the additional fields xxx and xxxlabel. See FT_DATATYPE_VOLUME

% for further details.

%

% Required fields:

% - dim, transform

%

% Optional fields:

% - brain, skull, scalp, gray, white, csf, or any other field with dimensions that are consistent with dim

% - unit, coordsys, fid

%

% Deprecated fields:

% - none

%

% Obsoleted fields:

% - none

%

% Revision history:

% (2012/latest) The explicit distunction between the indexed and probabilistic

% representation was made. For the indexed representation the additional

% xxxlabel cell-array was introduced.

%

% (2005) The initial version was defined.

%

% See also FT_DATATYPE, FT_DATATYPE_VOLUME, FT_DATATYPE_PARCELLATION

% Copyright (C) 2012, Robert Oostenveld

%

% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org

% for the documentation and details.

%

% FieldTrip is free software: you can redistribute it and/or modify

% it under the terms of the GNU General Public License as published by

% the Free Software Foundation, either version 3 of the License, or

% (at your option) any later version.

%

% FieldTrip is distributed in the hope that it will be useful,

% but WITHOUT ANY WARRANTY; without even the implied warranty of

% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the

% GNU General Public License for more details.

%

% You should have received a copy of the GNU General Public License

% along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.

%

% $Id$

% get the optional input arguments, which should be specified as key-value pairs

version = ft_getopt(varargin, 'version', 'latest');

segmentationstyle = ft_getopt(varargin, 'segmentationstyle'); % can be indexed or probabilistic

hasbrain = ft_getopt(varargin, 'hasbrain', 'no'); % no means that it is not required, if present it won't be removed

% convert from string into boolean

hasbrain = istrue(hasbrain);

if strcmp(version, 'latest')

segversion = '2012';

volversion = 'latest';

clear version

else

segversion = version;

volversion = version;

clear version

end

if isempty(segmentation)

return;

end

switch segversion

case '2012'

% convert the inside/outside fields, they should be logical rather than an index

if isfield(segmentation, 'inside')

segmentation = fixinside(segmentation, 'logical');

end

% make a list of fields that possibly represent a segmentation

fn = fieldnames(segmentation);

fn = setdiff(fn, 'inside'); % exclude the inside field from any conversions

sel = false(size(fn));

for i=1:numel(fn)

sel(i) = (isnumeric(segmentation.(fn{i})) || islogical(segmentation.(fn{i}))) && numel(segmentation.(fn{i}))==prod(segmentation.dim(1:3));

end

% only consider numeric fields of the correct size

fn = fn(sel);

% determine whether the style of the input fields is probabilistic or indexed

[indexed, probabilistic] = determine_segmentationstyle(segmentation, fn, segmentation.dim(1:3));

% ignore the fields that do not contain a segmentation

sel = indexed | probabilistic;

fn = fn(sel);

indexed = indexed(sel);

probabilistic = probabilistic(sel);

% convert from an exclusive to cumulative representation

% this is only only for demonstration purposes

% for i=1:length(sel)

% segmentation.(fn{sel(i)}) = volumefillholes(segmentation.(fn{sel(i)}));

% end

[dum, i] = intersect(fn, {'scalp', 'skull', 'brain'});

if numel(i)==3

% put them in the preferred order

fn(i) = {'scalp', 'skull', 'brain'};

end

[dum, i] = intersect(fn, {'skin', 'skull', 'brain'}); % this is not likely

if numel(i)==3

% put them in the preferred order

fn(i) = {'skin', 'skull', 'brain'};

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% ensure that the segmentation is internally consistent

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if any(probabilistic)

segmentation = fixsegmentation(segmentation, fn(probabilistic), 'probabilistic');

end

if any(indexed)

segmentation = fixsegmentation(segmentation, fn(indexed), 'indexed');

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% convert the segmentation to the desired style

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if isempty(segmentationstyle)

% keep it as it is

elseif strcmp(segmentationstyle, 'indexed') && any(probabilistic)

segmentation = convert_segmentationstyle(segmentation, fn(probabilistic), segmentation.dim, 'indexed');

indexed(probabilistic) = true; % these are now indexed

probabilistic(probabilistic) = false; % these are now indexed

elseif strcmp(segmentationstyle, 'probabilistic') && any(indexed)

segmentation = convert_segmentationstyle(segmentation, fn(indexed), segmentation.dim, 'probabilistic');

probabilistic(indexed) = true; % these are now probabilistic

indexed(indexed) = false; % these are now probabilistic

end % converting between probabilistic and indexed

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% add the brain if requested

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if hasbrain

if all(indexed)

fn = fieldnames(segmentation);

sel = false(size(fn));

for i=1:numel(fn)

sel(i) = any(strcmp(fn, [fn{i} 'label']));

end

fn = fn(sel);

if numel(fn)>1

ft_error('cannot construct a brain mask on the fly; this requires a single indexed representation');

else

tissue = segmentation.(fn{1});

tissuelabel = segmentation.([fn{1} 'label']);

if ~any(strcmp(tissuelabel, 'brain'))

threshold = 0.5;

smooth = 5;

% ensure that the segmentation contains the brain mask, if not then construct it from gray+white+csf

if ~all(ismember({'gray' 'white' 'csf'}, tissuelabel))

ft_error('cannot construct a brain mask on the fly; this requires gray, white and csf');

end

gray = tissue==find(strcmp(tissuelabel, 'gray'));

white = tissue==find(strcmp(tissuelabel, 'white'));

csf = tissue==find(strcmp(tissuelabel, 'csf'));

brain = gray + white + csf;

clear gray white csf seg

brain = volumesmooth(brain, smooth, 'brain');

brain = volumethreshold(brain, threshold, 'brain');

% store it in the output

segmentation.brain = brain;

end % try to construct the brain

end % if single indexed representation

elseif all(probabilistic)

if ~isfield(segmentation, 'brain')

if ~all(isfield(segmentation, {'gray' 'white' 'csf'}))

ft_error('cannot construct a brain mask on the fly; this requires gray, white and csf');

end

threshold = 0.5;

smooth = 5;

% ensure that the segmentation contains the brain mask, if not then construct it from gray+white+csf tissue probability maps

fprintf('creating brainmask ... using the sum of gray, white and csf tpms\n');

brain = segmentation.gray + segmentation.white + segmentation.csf;

brain = volumesmooth(brain, smooth, 'brain');

brain = volumethreshold(brain, threshold, 'brain');

% store it in the output

segmentation.brain = brain;

end

else

ft_error('cannot construct a brain mask on the fly; this requires a uniquely indexed or probabilitic representation');

end

end % if hasbrain

case '2005'

% the only difference is that the indexed representation for xxx did not have the xxxlabel field prior to the 2012 version

fn = fieldnames(segmentation);

sel = ~cellfun(@isempty, regexp(fn, 'label$'));

segmentation = rmfield(segmentation, fn(sel));

% furthermore it corresponds to the oldest version of the volume representation

volversion = '2003';

otherwise

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ft_error('unsupported version "%s" for segmentation datatype', segversion);

end

% the segmentation is a special type of volume structure, so ensure that it also fulfills the requirements for that

segmentation = ft_datatype_volume(segmentation, 'version', volversion);

% For the pass through ft_datatype_volume it is perhaps necessary to remove

% the fields that are specific for the segmentation and add them later again.

% At this moment ft_datatype_volume nicely passes all fields, so there is no

% special handling of the segmentation fields needed.