% Version 1.0.0
% Date: 2023-11-27
function result = ILMF(swir, k, IsAC)
[row, col, ~] = size(swir);
result = zeros(row, col);
IsEmission = zeros(row, col);
swir(swir==0)=1e-10;
for i_column = 1 : col  %第b列
    %% 取第i列元素
    %     fprintf("第%d列\n", i_column);
    x_observe = swir(:,i_column,:);
    x_observe = permute(x_observe,[1 3 2]);
    %% 初始化参数列表
    i = 0;
    index = false(row, 1); % 初始化
    alpha = zeros(row, 1); % 初始化
    last_threshold = -9999;
    threshold = 9999;
    while( abs(threshold - last_threshold) > 1e-4 && i < 3 )
        i = i + 1;
        % fprintf("第%d次迭代\n", i);
        %% 计算albedo factor
        junzhi = mean(x_observe, 'omitnan')';
        if(IsAC)
            albedo_factor = (x_observe *junzhi) ./ (junzhi'*junzhi);
            albedo_factor(index,:) = ( x_observe(index,:)  ./ exp(- k * alpha(index,:)')' ) * junzhi ./ (junzhi'*junzhi);
        end
        x_background = x_observe .* exp(- k * alpha')'; % 背景像素值
        % x_background = x_observe;
        % x_background(alpha > threshold,:) = []; % 把潜在的增强值像素剔除，用来计算更准的μ和C
        %% 计算分子分母
        y_n = log(x_background);
        % mu = mean(y_n, 'omitnan')';
        mu = mytrimmean(y_n, 5, 5)';
        Cb = cov(y_n, 'omitrows')';
        inv_Cb = pinv(Cb)';
        if(IsAC)
            ratio = mu - log(x_observe ./ repmat(albedo_factor, 1, length(k)))';
        else
            ratio = mu - log(x_observe)';
        end
        fenzi = ((k' * inv_Cb * ratio))'; %  用 mu 来估算 log(xr)
        fenmu = k' * inv_Cb * k;
        alpha = fenzi ./ fenmu;
        %% 计算阈值
        last_threshold = threshold;
        threshold = 2 * std(alpha, 'omitnan');
        index = (alpha > threshold);
        %         disp(threshold)
        % last_threshold = threshold;
    end
    alpha(alpha < 0 ) = 0;
    IsEmission(:, i_column) = index;
    result(:, i_column) = alpha;
    %     if(i == 30)
    %         disp(i_column)
    %     end
end
result = reshape(result, row, col);
% IsEmission = reshape(IsEmission, row, col);
%% 绘图
% figure('Position', [300 300 800 400])
% subplot(121)
% image(result)
% axis square
% subplot(122)
% imagesc(IsEmission)
% axis square
end