function [source] = ft_datatype_source(source, varargin)

% FT_DATATYPE_SOURCE describes the FieldTrip MATLAB structure for data that is

% represented at the source level. This is typically obtained with a beamformer of

% minimum-norm source reconstruction using FT_SOURCEANALYSIS.

%

% An example of a source structure obtained after performing DICS (a frequency domain

% beamformer scan) is shown here

%

% pos: [6732x3 double] positions at which the source activity could have been estimated

% inside: [6732x1 logical] boolean vector that indicates at which positions the source activity was estimated

% dim: [xdim ydim zdim] if the positions can be described as a 3D regular grid, this contains the

% dimensionality of the 3D volume

% cumtapcnt: [120x1 double] information about the number of tapers per original trial

% time: 0.100 the latency at which the activity is estimated (in seconds)

% freq: 30 the frequency at which the activity is estimated (in Hz)

% pow: [6732x120 double] the estimated power at each source position

% powdimord: 'pos_rpt' defines how the numeric data has to be interpreted,

% in this case 6732 dipole positions x 120 repetitions (i.e. trials)

% cfg: [1x1 struct] the configuration used by the function that generated this data structure

%

% Required fields:

% - pos

%

% Optional fields:

% - inside, pow, coh, eta, mom, ori, leadfield, filter, or any other field with dimensions that are consistent with pos or dim

% - dim, transform, unit, coordsys, time, freq, cumtapcnt, dimord

%

% Deprecated fields:

% - method, outside

%

% Obsoleted fields:

% - xgrid, ygrid, zgrid, transform, latency, frequency

%

% Revision history:

%

% (2014) The subfields in the avg and trial fields are now present in the

% main structure, e.g. source.avg.pow is now source.pow. Furthermore, the

% inside is always represented as logical vector.

%

% (2011) The source representation should always be irregular, i.e. not

% a 3-D volume, contain a "pos" field and not contain a "transform".

%

% (2010) The source structure should contain a general "dimord" or specific

% dimords for each of the fields. The source reconstruction in the avg and

% trial substructures has been moved to the toplevel.

%

% (2007) The xgrid/ygrid/zgrid fields have been removed, because they are

% redundant.

%

% (2003) The initial version was defined

%

% See also FT_DATATYPE, FT_DATATYPE_COMP, FT_DATATYPE_DIP, FT_DATATYPE_FREQ,

% FT_DATATYPE_MVAR, FT_DATATYPE_RAW, FT_DATATYPE_SOURCE, FT_DATATYPE_SPIKE,

% FT_DATATYPE_TIMELOCK, FT_DATATYPE_VOLUME

% Copyright (C) 2013-2014, 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$

% FIXME: I am not sure whether the removal of the xgrid/ygrid/zgrid fields

% was really in 2007

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

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

if strcmp(version, 'latest') || strcmp(version, 'upcoming')

version = '2014';

end

if isempty(source)

return;

end

% old data structures may use latency/frequency instead of time/freq. It is

% unclear when these were introduced and removed again, but they were never

% used by any FieldTrip function itself.

if isfield(source, 'frequency')

source.freq = source.frequency;

source = rmfield(source, 'frequency');

end

if isfield(source, 'latency')

source.time = source.latency;

source = rmfield(source, 'latency');

end

switch version

case '2014'

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

% ensure that it has individual source positions

source = fixpos(source);

if isfield(source, 'inside')

% ensure that it is always logical

source = fixinside(source, 'logical');

end

if isfield(source, 'coordsys')

% ensure that it is in lower case

source.coordsys = lower(source.coordsys);

end

if isfield(source, 'unit')

% ensure that it is in lower case

source.unit = lower(source.unit);

end

% remove obsolete fields

if isfield(source, 'method')

source = rmfield(source, 'method');

end

if isfield(source, 'transform')

source = rmfield(source, 'transform');

end

if isfield(source, 'xgrid')

source = rmfield(source, 'xgrid');

end

if isfield(source, 'ygrid')

source = rmfield(source, 'ygrid');

end

if isfield(source, 'zgrid')

source = rmfield(source, 'zgrid');

end

if isfield(source, 'avg') && isstruct(source.avg) && isfield(source, 'trial') && isstruct(source.trial) && ~isempty(intersect(fieldnames(source.avg), fieldnames(source.trial)))

% it is not possible to convert both since they have the same field names

ft_warning('removing ''avg'', keeping ''trial''');

source = rmfield(source, 'avg');

end

if isfield(source, 'avg') && isstruct(source.avg)

% move the average fields to the main structure

fn = fieldnames(source.avg);

for i=1:length(fn)

dat = source.avg.(fn{i});

if isequal(size(dat), [1 size(source.pos,1)])

source.(fn{i}) = dat';

else

source.(fn{i}) = dat;

end

clear dat

end % j

source = rmfield(source, 'avg');

end

if isfield(source, 'trial') && isstruct(source.trial)

npos = size(source.pos,1);

% concatenate the fields for each trial and move them to the main structure

fn = fieldnames(source.trial);

for i=1:length(fn)

% some fields are descriptive and hence identical over trials

if any(strcmp(fn{i}, {'csdlabel' 'label' 'filterdimord' 'leadfielddimord'}))

source.(fn{i}) = source.trial(1).(fn{i});

continue

end

% start with the first trial

dat = source.trial(1).(fn{i});

datsiz = getdimsiz(source, fn{i});

nrpt = datsiz(1);

datsiz = datsiz(2:end);

if iscell(dat)

datsiz(1) = nrpt; % swap the size of pos with the size of rpt

val = cell(npos,1);

if isfield(source, 'inside')

indx = find(source.inside);

else

indx = 1:npos;

end

for k=1:length(indx)

val{indx(k)} = nan(datsiz);

val{indx(k)}(1,:,:,:) = dat{indx(k)};

end

% concatenate all data as {pos}_rpt_etc

for j=2:nrpt

dat = source.trial(j).(fn{i});

for k=1:length(indx)

val{indx(k)}(j,:,:,:) = dat{indx(k)};

end

end % for all trials

source.(fn{i}) = val;

else

% concatenate all data as pos_rpt_etc

val = nan([datsiz(1) nrpt datsiz(2:end)]);

val(:,1,:,:,:) = dat(:,:,:,:);

for j=2:length(source.trial)

dat = source.trial(j).(fn{i});

val(:,j,:,:,:) = dat(:,:,:,:);

end % for all trials

source.(fn{i}) = val;

% else

% siz = size(dat);

% if prod(siz)==npos

% siz = [npos nrpt];

% elseif siz(1)==npos

% siz = [npos nrpt siz(2:end)];

% end

% val = nan(siz);

% % concatenate all data as pos_rpt_etc

% val(:,1,:,:,:) = dat(:);

% for j=2:length(source.trial)

% dat = source.trial(j).(fn{i});

% val(:,j,:,:,:) = dat(:);

% end % for all trials

% source.(fn{i}) = val;

end

end % for each field

source = rmfield(source, 'trial');

end % if trial

% ensure that it has a dimord (or multiple for the different fields)

source = fixdimord(source);

if isfield(source, 'inside')

% ensure that for positions outside the brain it is [], not nan

if isfield(source, 'leadfield')

source.leadfield(~source.inside) = {[]};

end

if isfield(source, 'filter')

source.filter(~source.inside) = {[]};

end

end

if isfield(source, 'leadfield') && ~isfield(source, 'label') && isfield(source, 'cfg')

% try to determine the channel labels from the cfg

label = ft_findcfg(source.cfg, 'channel');

if ~isempty(label)

source.label = label;

end

end

if isfield(source, 'filter') && ~isfield(source, 'label') && isfield(source, 'cfg')

% try to determine the channel labels from the cfg

label = ft_findcfg(source.cfg, 'channel');

if ~isempty(label)

source.label = label;

end

end

% ensure that all data fields have the correct dimensions

fn = getdatfield(source);

for i=1:numel(fn)

dimord = getdimord(source, fn{i});

dimtok = tokenize(dimord, '_');

dimsiz = getdimsiz(source, fn{i}, numel(dimtok));

if numel(dimsiz)>=3 && strcmp(dimtok{1}, 'dim1') && strcmp(dimtok{2}, 'dim2') && strcmp(dimtok{3}, 'dim3')

% convert it from voxel-based representation to position-based representation

try

source.(fn{i}) = reshape(source.(fn{i}), [prod(dimsiz(1:3)) dimsiz(4:end) 1]);

catch

ft_warning('could not reshape %s to the expected dimensions', fn{i});

end

end

end

% ensure that the structure has all required fields

for required={'pos'}

assert(isfield(source, required), 'required field "%s" is missing', required{:});

end

case '2011'

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

% ensure that it has individual source positions

source = fixpos(source);

% remove obsolete fields

if isfield(source, 'xgrid')

source = rmfield(source, 'xgrid');

end

if isfield(source, 'ygrid')

source = rmfield(source, 'ygrid');

end

if isfield(source, 'zgrid')

source = rmfield(source, 'zgrid');

end

if isfield(source, 'transform')

source = rmfield(source, 'transform');

end

% ensure that it has a dimord (or multiple for the different fields)

source = fixdimord(source);

case '2010'

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

% ensure that it has individual source positions

source = fixpos(source);

% remove obsolete fields

if isfield(source, 'xgrid')

source = rmfield(source, 'xgrid');

end

if isfield(source, 'ygrid')

source = rmfield(source, 'ygrid');

end

if isfield(source, 'zgrid')

source = rmfield(source, 'zgrid');

end

% ensure that it has a dimord (or multiple for the different fields)

source = fixdimord(source);

case '2007'

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

% ensure that it has individual source positions

source = fixpos(source);

% remove obsolete fields

if isfield(source, 'dimord')

source = rmfield(source, 'dimord');

end

if isfield(source, 'xgrid')

source = rmfield(source, 'xgrid');

end

if isfield(source, 'ygrid')

source = rmfield(source, 'ygrid');

end

if isfield(source, 'zgrid')

source = rmfield(source, 'zgrid');

end

case '2003'

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

if isfield(source, 'dimord')

source = rmfield(source, 'dimord');

end

if ~isfield(source, 'xgrid') || ~isfield(source, 'ygrid') || ~isfield(source, 'zgrid')

if isfield(source, 'dim')

minx = min(source.pos(:,1));

maxx = max(source.pos(:,1));

miny = min(source.pos(:,2));

maxy = max(source.pos(:,2));

minz = min(source.pos(:,3));

maxz = max(source.pos(:,3));

source.xgrid = linspace(minx, maxx, source.dim(1));

source.ygrid = linspace(miny, maxy, source.dim(2));

source.zgrid = linspace(minz, maxz, source.dim(3));

end

end

otherwise

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

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

end

function pos = grid2pos(xgrid, ygrid, zgrid)

[X, Y, Z] = ndgrid(xgrid, ygrid, zgrid);

pos = [X(:) Y(:) Z(:)];

function pos = dim2pos(dim, transform)

[X, Y, Z] = ndgrid(1:dim(1), 1:dim(2), 1:dim(3));

pos = [X(:) Y(:) Z(:)];

pos = ft_warp_apply(transform, pos, 'homogenous');