%% 读取云掩膜结果
function Cloud = DetectCloud(fpath, flist, xmlStruct)
Cloud_file = fullfile(fpath, 'Cloud.mat');
if(exist(Cloud_file, 'file'))
    load(Cloud_file, 'Cloud')
else
    %% 判断卫星是高分系列还是资源系列
    series = xmlStruct.SatelliteID(1:2); % 卫星系列（GF/ZY）
    if(strcmp(series, 'GF'))
        index = {11:20; 30:60; 192; 270:272};
        load('data/mat_data_GF5B/Esun.mat', 'Esun');
    elseif(strcmp(series, 'ZY'))
        index = {5:10; 15:30; 96; 135:136};
        load('data/mat_data_ZY1E/Esun.mat', 'Esun');
    else
        fprintf("请检查xmlStruct.SatelliteID是否准确\n");
    end

    %% 读取所需数据
    mydisp('读取SWIR数据...')
    swir = readgeoraster(flist.SWIR);
    sw_rad_cal = Read_Raw(flist.swir_rad, 1); % 辐射定标系数 radiance calculation
    swir = Rad_Calibration(swir, sw_rad_cal); % 短波红外辐射定标  单位：W·m-2·sr-1·μm-1
    mydisp('读取VNIR数据...')
    vnir = readgeoraster(flist.VNIR);
    vn_rad_cal = Read_Raw(flist.vnir_rad, 1); % 辐射定标系数
    vnir = Rad_Calibration(vnir, vn_rad_cal); % 可见光与近红外辐射定标
    data = MergeData(vnir, swir);
    FWHM_SWIR = Read_Raw(flist.swir_fwhm, 2); % FWHM 包含central wavelength 中心波长 和 full width at half maxima 半宽
    FWHM_VNIR = Read_Raw(flist.vnir_fwhm, 2);
    fwhm = cat(1, FWHM_VNIR, FWHM_SWIR);

    %% 数据准备
    date = datetime(xmlStruct.StartTime);
    Day1 = juliandate(date);
    Day2 = juliandate(datetime(year(date), 1, 1));
    Days = Day1 - Day2;
    D = (1 + 0.0167 * sin(2 * pi * (Days - 93.5) / 365)); % 日地平均距离（天文单位）
    Esun = reshape(Esun, 1, 1, []);
    Esun = repmat(Esun, size(swir, 1) , size(swir, 2));
    SolarZenith = str2double(xmlStruct.SolarZenith); % 太阳天顶角

    %% 合并为对应的一个宽波段数据
    % T1波段 433~471nm (第11~第20波段)
    T1 = mergeband(data, fwhm, index{1}, D, Esun, SolarZenith);

    % T2波段 510~640nm (第30~第60波段)
    T2 = mergeband(data, fwhm, index{2}, D, Esun, SolarZenith);

    % T3波段 1350~1360nm (第192波段)   卷云波段(1360~1390 nm)
    T3 = mergeband(data, fwhm, index{3}, D, Esun, SolarZenith);

    % T4波段 2000nm (第270~第272波段)
    T4 = mergeband(data, fwhm, index{4}, D, Esun, SolarZenith);
    %% 计算云掩膜并编码
    c1 = double((T1 > 0.3) & (T2 > 0.3)); % 厚云
    c2 = double((T3 > 0.2) & (T1 > 0.15)); % 薄云
    c3 = double((T1 ./ T4 > 10) & (T4 ./ T3 < 1)); % 云影
    Cloud = uint8(c1 * 4 + c2 * 2 + c3);
    %% 后处理
    Cloud = Cloud_post(Cloud, 15);
    imagesc(Cloud);
    %%
    imagesc(Cloud)
    save(Cloud_file, 'Cloud');
    mydisp('已保存Cloud文件')
end
end
%% 函数
% 合并为宽波段数据
function T = mergeband(data, fwhm, index, D, Esun, SolarZenith)
fwhm = fwhm(index,2);
rho = (pi * data(:, :, index) * D^2) ./ (Esun(:, :, index) * cos(pi * SolarZenith / 180)); % rho为大气顶层的表观反射率(各像素，各波段)
w = reshape(fwhm, 1, 1, [])/sum(fwhm);
w = repmat(w, size(data, 1), size(data, 2));
T = sum(rho .* w, 3);
end
%% 绘图
% figure
% subplot(2, 4, 1)
% imagesc(T1)
% colorbar
% axis square
% subplot(2, 4, 2)
% imagesc(T2)
% colorbar
% axis square
% subplot(2, 4, 3)
% imagesc(T3)
% colorbar
% axis square
% subplot(2, 4, 4)
% imagesc(T4)
% colorbar
% axis square
% subplot(2, 4, 5)
% imagesc(T1>0.3)
% colorbar
% axis square
% subplot(2, 4, 6)
% imagesc(T2>0.3)
% colorbar
% axis square
% subplot(2, 4, 7)
% imagesc((T3 > 0.04) & (T1 > 0.15))
% colorbar
% axis square
% subplot(2, 4, 8)
% imagesc((T1 ./ T4 > 7.5) & (T4 ./ T3 < 1))
% colorbar
% axis square
