function [varargout] = barumerli2021_metrics(varargin)
%barumerli2021_metrics - extract localization metrics
%   Usage: [mean_error, bias] = barumerli2021_metrics(doa, parameters) 
%
%   Input parameters:
%
%   doa: Struct in returned from the barumerli2021's model with
%        estimated and real directions of arrival
%
%   BARUMERLI2021_METRICS(...) returns psychoacoustic performance 
%   parameters for experimental response patterns. 
%   doa is a struct where actual and estimated directions of arrival must
%   be provided. If no input params are provided the returned metrics
%   resemble the ones provided in the original paper, see Reijiners et al. (2014).
%   This script is a wrapper for localizationerror.
%
%   If parameter is provided, barumerli2021_metric is a wrapper for localizationerror
%   with the parameter as the localization error. 
%
%   Url: http://amtoolbox.org/amt-1.1.0/doc/modelstages/barumerli2021_metrics.php

% Copyright (C) 2009-2021 Piotr Majdak, Clara Hollomey, and the AMT team.
% This file is part of Auditory Modeling Toolbox (AMT) version 1.1.0
%
% This program 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.
%
% This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.

%   See also: demo_barumerli2021 barumerli2021 localizationerror
%
%   References: katz2014


%   Author: Roberto Barumerli
%   Dept. of Information Engineering, University of Padova

% assert(isa(doa.coords, 'barumerli2021_coordinates') && isa(doa_real, 'barumerli2021_coordinates'),  'barumerli2021_metrics: parameters should be wrapped into the barumerli2021_coordinates class')
% assert(size(doa.estimations, 1) == doa_real.count_positions, 'barumerli2021_metrics: size mismatch')

% parameters
% empty                     plot stuff like reijniers with mean spherical error
% m, 'middle_metrics'       return middlebrooks metrics as a struct
% doa, doa_real, 'm'        return m matrix
% doa, doa_real, <error>    compute error from localizationerror.m
% m, <error>                as above

if strcmp(varargin{2}, 'middle_metrics')
    assert(size(varargin{1}, 2) == 8, 'Please provide m matrix')
    m = varargin{1};
    % lateral_bias 
    exp.accL = localizationerror(m, 'accL'); 
    % lateral_rms_error
    exp.rmsL = localizationerror(m, 'rmsL'); 
    % elevation_bias
    exp.accP = localizationerror(m, 'accP'); 
    % local_rms_polar
    exp.rmsP = localizationerror(m, 'rmsPmedianlocal'); 
    % quadrant_err
    exp.querr = localizationerror(m, 'querrMiddlebrooks'); 
    varargout{1} = exp;
elseif strcmp(varargin{3}, 'm')
    assert(isfield(varargin{1}, 'estimations') & isa(varargin{2}, 'barumerli2021_coordinates'), 'If looking for m matrix please give doa as a struct and doa_real as barumerli2021_coordinates object(see barumerli2021)')
    varargout{1} = local_returnmatrixlocalizationerror(varargin{1}, varargin{2});
else
    if isfield(varargin{1}, 'estimations') && isa(varargin{2}, 'barumerli2021_coordinates')
        m = local_returnmatrixlocalizationerror(varargin{1}, varargin{2});
        errorflag = varargin{3};
    elseif size(varargin{1}, 2) == 8
        m = varargin{1};
        errorflag = varargin{2};
    else
        error('something went wrong!')
    end

    [varargout{1}, meta, par] = localizationerror(m, errorflag);
    
    if length(varargout) > 1
        varargout{2}=meta;
    end
    if length(varargout) > 2
        varargout{3}=par; 
    end

end 

function m = local_returnmatrixlocalizationerror(doa, doa_real)
    assert(size(doa.estimations, 3) == 3)
   
    doa_est_cart = barumerli2021_coordinates(reshape(doa.estimations, [], 3), 'cartesian');
    
    %% compute the metric relying on localizationerror.m
    doa_real_sph = doa_real.return_positions('spherical');
    doa_est_sph = doa_est_cart.return_positions('spherical');
    doa_real_hor = doa_real.return_positions('horizontal-polar');
    doa_est_hor = doa_est_cart.return_positions('horizontal-polar');
    
    num_rep = size(doa_est_cart.pos, 1)/size(doa_real.pos, 1);
    
    m = zeros(size(doa_real.pos, 1)*num_rep, 8);
    m(:, 1:2) = repmat(doa_real_sph(:, [1 2]), num_rep, 1);
    m(:, 3:4) = doa_est_sph(:, [1 2]);
    m(:, 5:6) = repmat(doa_real_hor(:,[1 2]), num_rep, 1);
    m(:, 7:8) = doa_est_hor(:, [1 2]);