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

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

%%
function alpha = Filter(L, s, Niter, IsAC)
%% 输出参数
% L：行数 * 波段数
% s：波段数 * 1
% Niter: 迭代次数
param1 = 5;
param2 = 5;
%%
mu = cell(Niter, 1);
C = cell(Niter, 1);
alpha = cell(Niter, 1);
w = cell(Niter, 1);
%% 初始化
N = length(L); % 行数
% mu{1} = trimmean(L, param1)'; % 波段数 * 1
mu{1} = mytrimmean(L, param1, param2)';
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)); % 超过平均不计算入背景光谱
    %     nan_idx = isnan(d);
    %     nan_idx = logical(sum(nan_idx, 2));
    %     d(emission_idx,:) = [];
    mu{k} = mytrimmean(d, param1, param2)';
    tempA = mu{k} .* s; % 8 * 1
    LC = L - r .* alpha{k - 1} * tempA'; % 2056 * 8

    %%
    %     % 计算每行数据的平均值
    %     mean_LC = mean(LC, 2);
    %     % 对平均值进行排序
    %     [~, sorted_idx] = sort(mean_LC);
    %     % 计算要排除的像素的数量
    %     num_rows = size(LC, 1);
    %     num_rows_to_exclude = floor(0.05 * num_rows);
    %     % 计算需要保留的行的索引范围
    %     keep_indices = (num_rows_to_exclude+1):(num_rows - num_rows_to_exclude);
    %     % 保留 LC 中符合条件的行
    %     LC_cleaned = LC(sorted_idx(keep_indices), :);
    %     % 计算协方差矩阵
    %     C{k} = cov(LC_cleaned, 'omitrows');
    C{k} = cov(LC, 'omitrows');
    %%
    %     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