function [X,it] = hyperbaseline(Y,s,alpha,beta,tol,pcgitmax)
%Estimation de la ligne de base d'une image hyperspectrale.
%
% Utilisation:
%
%   [X, it] = robuste3(Y, s, alpha, beta, tol, pcgitmax)
%
% Entrées:
%   Y = image hyperspecrale dont les 2 premières dimensions sont spatiales
%       et la 3ème dimension est spectrale
%   s = paramètre seuil de la fonction de Huber
%   alpha = paramètre de régularisation spectrale
%   beta  = paramètre de la régularisation spatiale
%   [tol] = tolérance sur l'erreur relative pour convergence (1e-6)
%   [pcgitmax] = nombre max d'itération de l'algorithme PCG (1000)
%
% Sorties:
%   X = signal sans ligne de base
%   it = nombre d'itérations pour converger.
if nargin<5
    tol = 1e-6;
end
if nargin<6
    pcgitmax = 1000;
end

[m,n,p] = size(Y);
y = Y(:);

% les operateurs de dérivée
Dv = spdiags([ones(m,1) -ones(m,1)],[0 1],m-1,m);   % dim 1 = spatiale
Dh = spdiags([ones(n,1) -ones(n,1)],[0 1],n-1,n);   % dim 2 = spatiale
Dl = spdiags([ones(p,1) -ones(p,1)],[0 1],p-1,p);   % dim 3 = spectrale

% matrice à inverser
A = speye(N)/(2*c) + alpha*kron(Dl'*Dl,speye(n*m)) + ...
    beta*kron(speye(n*p),Dv'*Dv) + ...
    beta*kron(speye(p),kron(Dh'*Dh,speye(m)));

% minimisation semi-quadratique
c = 0.99 * 1/2;

% début
it = 0;             % compteur d'itérations
x1 = y;             % estimations précédentes
xp=ones(size(y));

while sum((x1-xp).^2)/sum(xp.^2) > tol,
    
    it = it + 1;
    xp = x1;
    res = y - x1;       %  résidus
    
    % estimer z
    z = (res*(2*c-1)) .* (res<s) + (c*2*s-res) .* (res>=s);
    
    % estimer x
    x = pcg(A,y+z,tol,pcgitmax);

    x1 = x;
end

X = reshape(x,[m n p]);
