% Version: 1.1.0
% Reviser: Huayiwang
% Feature: Robust Matched Filter
% Date: 2023-12-02

function result = RMF_all(swir, s, IsAC, Trim_pct)
% IsAC: albedo factor
[r,c,w] = size(swir);
Niter = 2; % 迭代次数
result = {zeros(size(swir, [1 2]))};
result = repmat(result, Niter, 1);
L = reshape(swir,[r*c,w]); % 取出一列元素
alpha = Filter(L, s, Niter, IsAC, Trim_pct);
alpha = alpha{Niter};
result = reshape(alpha, [r, c]);
end

%%
function alpha = Filter(L, s, Niter, IsAC, Trim_pct)
%% 输出参数
% 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 = pinv(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;
%%
IsWeight = false; % True：需要L1权重，反演结果干净，偏小；false：不需要L1权重，反演结果较杂乱，难裁剪
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
    d = L - r .* alpha{k - 1} * mu{k - 1}' .* s';
    id_temp = alpha{k - 1} > 0;
    emission_idx = alpha{k - 1} > mean(alpha{k - 1}(id_temp)); % 超过平均不计算入背景光谱
    d(emission_idx,:) = [];
    mu{k} = trimmean(d, Trim_pct)';
    tempA = mu{k} .* s; % 8 * 1
    LC = L - r .* alpha{k - 1} * tempA'; % 2056 * 8
    C{k} = cov(LC, 'omitrows'); % 8 * 8
    tempC = pinv(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