% Version: 1.0.0
% Reviser: Huayiwang
% Feature: Match according to the column

function [result, likehold] = MF(swir, k)

[row, col, ~] = size(swir);
result = zeros(row, col);
likehold = zeros(row, col);
% 进度条输出
% interval = 20;
% percent_out = -1 * interval; % 初始化
for i = 1 : col %第b列   
%     percent = round(i / col * 100);
%     if(percent == percent_out + interval)
%         percent_out = percent;
%         fprintf('%s 反演进度为：%d%%\n', datetime('now'), percent_out);
%     end
    xt_obs = swir(:,i,:); % 取出一列元素
    xt_obs = permute(xt_obs,[1 3 2]); % 去除列向量层级，变为行*波段
    
    % 真实的数据中，可能会出现 某个波段所有像素 全是0 的情况
    index_all_zero = find(all(xt_obs == 0, 1));% 看看第几个波段 全是0
    k_temp = k;
    k_temp(index_all_zero,:) = [];
    xt_obs(:,index_all_zero) = [];
    % 真实的数据中，可能会出现 个别像素个别波段 是0 的情况，求对数后影响求逆
    index_zero = any(xt_obs==0, 2); % 看看第几个像素 含0
    xt_obs(index_zero, :) = nan;
    
    if sum(all(isnan(xt_obs))) == 0 % 如果temp全为空，就不计算
        %25:38 2230-2330nm      改为  8:55 2110-2448nm
        mean_data = nanmean(xt_obs)';
        albedo_factor = (xt_obs * mean_data) ./ (mean_data' * mean_data);
        mu = nanmean(xt_obs)';
        Cb = nancov(xt_obs)';
        Cb_inv = pinv(Cb);
        fenzi = (mu - xt_obs')' * Cb_inv * (mu .* k_temp);
        fenmu = albedo_factor .* (mu .* k_temp)' * Cb_inv * (mu .* k_temp);
        alpha = fenzi ./ fenmu;
        alpha(alpha < 0 ) = 0;
%         mut = mu .* exp(-k_temp * alpha'); % mu - k_temp * alpha';
%         Ct = nancov(current_col_swir .* exp(-k_temp * alpha')' );
%         Ct(isnan(Ct)) = 0;
%         Ct(isinf(Ct)) = 0;
%         Ct_inv = pinv(Ct);
%         
%         Pt = exp(-1/2 .* ( (mut - current_col_swir')' * Ct_inv * (mut - current_col_swir')) ) ./ (  sqrt(det(Ct)));
%         Pb = exp(-1/2 .* ( (mu - current_col_swir')' * Cb_inv * (mu - current_col_swir')) ) ./ (  sqrt(det(Cb)));
%         likehold(:,i) = diag(Pt ./ Pb); % 越大越好
        result(:,i) = alpha;
    end
end
result(isnan(result) & isinf(result)) = 0;
result = reshape(result,row,col);
end