% Version: 1.0
% Reviser: Huayiwang
% Feature: GLRT Matched Filter
% Date: 2024-4-2

function result = GLRTMF(swir, s)
Niter = 1; % 迭代次数
result = {zeros(size(swir, [1 2]))};
result = repmat(result, Niter, 1);
for j = 1 : size(swir, 2) % 按列进行处理
    L = swir(:, j, :); % 取出一列元素
    L = squeeze(L); % 去除列向量层级，变为行*波段
    alpha = Filter(L, s, Niter);
    for i = 1 : Niter
        result{i}(:, j) = alpha{i}; % 取最后一次迭代结果
    end
end
result = result{1};
end

%%
function alpha = Filter(L, s, Niter)
%% 输出参数
% L：行数 * 波段数
% s：波段数 * 1
% Niter: 迭代次数
%%
mu = cell(Niter, 1);
C = cell(Niter, 1);
alpha = cell(Niter, 1);
%% 初始化
N = length(L); % 行数
mu{1} = mean(L, 'omitnan')'; % 波段数 * 1
C{1} = cov(L, 'omitrows'); % 波段数 * 波段数
alpha{1} = zeros(N, 1);
%% Albedo Factor
tempA = mu{1} .* s; % 8 * 1
tempC = C{1} \ tempA; % 8 * 1
tempB = mu{1} - L'; % 8 * 2056
alpha{1} = (tempB' * tempC + sum(s)) ./ (tempA' * tempC + (s .* tempA)' * C{1} * tempB)';
alpha{1}(alpha{1} < 0) = 0;

%%
% for k = 2 : Niter
%     mu{k} = mean(L - alpha{k - 1} * mu{k - 1}' .* s', 'omitnan')';
%     tempA = mu{k} .* s; % 8 * 1
%     LC = L - alpha{k - 1} * tempA'; % 2056 * 8
%     C{k} = pinv(cov(LC, 'omitrows')); % 8 * 8
%     tempC = C{k} * tempA; % 8 * 1
%     alpha{k} = zeros(N, 1);
%     tempB = mu{k} - L'; % 8 * 2056
%     alpha{k} = (tempB' * tempC + sum(s)) ./ (tempA' * tempC + (s .* tempA)' * C{k} * tempB)';
%     alpha{k}(alpha{k} < 0) = 0;
% end
end