% Version: 1.0
% Reviser: Huayiwang
% Feature: Iterative Logarithmic Matched Filter
% Date: 2024-4-02

function result = LMF(swir, s, IsWeight, IsAC)
% IsAC: albedo factor

Niter = 8; % 迭代次数
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, IsWeight, IsAC);
    for i = 1 : Niter
        result{i}(:, j) = alpha{i}; % 取最后一次迭代结果
    end
end
end

%%
function alpha = Filter(L, s, Niter, IsWeight, IsAC)
%% 输出参数
% L：行数 * 波段数
% s：波段数 * 1
% Niter: 迭代次数
%%
mu = cell(Niter, 1);
C = cell(Niter, 1);
alpha = cell(Niter, 1);
w = 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
if(IsAC)
    r = (L * mu{1}) / (mu{1}' * mu{1}); % 2056 * 1
else
    r = 1;
end
tempB = mu{1} - L'; % 8 * 2056
alpha{1} = tempB' * tempC ./ (r * tempA' * tempC);
alpha{1}(alpha{1} < 0) = 0;

%%
for k = 2 : Niter
    if(IsWeight) % 是否需要权重平衡
        epsilon = 1e-10;
        w{k} = 1 ./ (alpha{k - 1} * 2 + epsilon); % 2056 * 1
    else
        w{k} = zeros(N, 1);
    end
    mu{k} = mean(L - r .* alpha{k - 1} * mu{k - 1}' .* s', 'omitnan')';
    tempA = mu{k} .* s; % 8 * 1
    LC = L - r .* alpha{k - 1} * tempA'; % 2056 * 8
    C{k} = 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 - w{k}) ./ (r * tempA' * tempC);
    alpha{k}(alpha{k} < 0) = 0;
end
end