function [XN,factors,params]=NormaliseGNAT(X,ppmH,method,features,value,peakID)
% Description:
%       Normalize spectral matrix X to total area, single feature, or integral of set of features.
%
% Input:
%       X        : stack of 1D spectra
%       ppmH     : chemical shift vector
%       method   : 'total' for Total Area,
%                  'largest' for the height of the highest peaks,
%                  'PQN' for Probablistic Quotient,
%                  'quantile' for Quantile Normalization,
%                  'intensity' for normalization to single feature,
%                  'integral' for normalization to sum of set of features
%       features : only required for 'intensity' or 'integral'.
%                  For 'intensity', the ppm of the feature in X to normalize to - e.g. [10] for X(:,10).
%                  For 'integral', the range of features in X that span the peak to normalize to - eg [-.05,0.05]
%
%                  If method is set to 'quantile',
%                  optionally set feature to 'median' to take median of the ranked values instead of the mean.
%
% Output:
%       XN      : N x P matrix of normalized data
%       factors : N calculated normalization factors: XN = X / factors
%
% Log:
%             'factors' output for PQN now includes the initial 'total' normalization factors
%             'factors' output for quantile returns nan to avoid error
% Example run:
% [XN,factors]=normalize(X,ppmR,'PQN');
%
%
%   This is a part of the GNAT
%   Copyright  2024  <Mathias Nilsson>%
%   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 2 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, write to the Free Software Foundation, Inc.,
%   51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
%
%   Dr. Mathias Nilsson
%   School of Chemistry, University of Manchester,
%   Oxford Road, Manchester M13 9PL, UK
%   Telephone: +44 (0) 161 306 4465
%   Fax: +44 (0)161 275 4598
%
%   Hugo da Silva Rocha, PhD Student
%   School of Chemistry, University of Manchester,
%   hugo.rocha@postgrad.manchester.ac.uk
%   Telephone: +44 (0) 7871080481

XN=zeros(size(X));

if nargin<2
    error('incorrect number of input parameters');
end

[obs dim]=size(X);

switch lower(method)
    case 'total'
        factors=repmat(NaN,[1 obs]);
        for i=1:obs
            factors(i)=sum(X(i,:));
            XN(i,:)=value*(X(i,:)./factors(i));
        end
    case 'largest'
        factors=max(max(X));
        XN=X/factors;
    case 'pqn'
        [X,factorsT] = NormaliseGNAT(X,ppmH,'total',features,value);      
        X = X.*100;                                    
        X(0==X)=0.00000001;
        normRef=repmat(median(X),size(X,1),1);
        F=X./(normRef);
        for i=1:obs
            factors(i)=median(F(i,:));
            XN(i,:)=X(i,:)./factors(i);
        end
        factors = factors .* factorsT;                
    case 'intensity'
        factors=X(:,peakID);
        for i=1:obs
            XN(i,:)=X(i,:)./factors(i);
        end
    case 'integral'
        [~,feature(1)]=min(abs(ppmH-features(1)));
        [~,feature(2)]=min(abs(ppmH-features(2)));
        for i=1:obs
            factors(i)=trapz(X(i,feature(2):feature(1)));
            XN(i,:)=value*(X(i,:)./factors(i));
        end
    case 'quantile'
        if exist('features','var')
            XN = transpose(quantilenorm(X',peakID,'true'));
        else
            XN = transpose(quantilenorm(X'));
        end
        factors = nan;      
    case 'none'
        norm_factor = max(max(X));
        XN = X*norm_factor;
        factors = 0;
end
end