clc;
clear;
close all;
tic;
%%
% 设置工作目录
rootDir = 'E:\2024.12 山东东营结果';

% 获取所有子目录
subDirs = dir(fullfile(rootDir, 'CCUS*'));
subDirs = subDirs([subDirs.isdir] & ~ismember({subDirs.name}, {'.', '..'}));

% 初始化变量
allRasters = {}; % 存储所有栅格对象
minLat = Inf; maxLat = -Inf;
minLon = Inf; maxLon = -Inf;

foldernum = length(subDirs);

GasType = 'CO2';
% GasType = 'CH4';
% GasType = 'RGB';
if(strcmp(GasType, 'CO2'))
    cmax = 50;
elseif(strcmp(GasType, 'CH4'))
    cmax = 1;
end
imagenum = 0;

for i = 1:foldernum
    % 进入子目录
    subdirPath = fullfile(rootDir, subDirs(i).name);

    % 读取当前子目录下所有以 GasType 开头且以 '.tif' 结尾的文件
    tifFiles = dir(fullfile(subdirPath, [GasType, '*.tif']));

    % 使用 cellfun 过滤掉文件名中包含 'thresh' 的文件
    tifFiles = tifFiles(~cellfun(@(x) contains(x, 'thresh'), {tifFiles.name}));

    % 使用 cellfun 过滤掉文件名中包含 'thresh' 的文件
    tifFiles = tifFiles(~cellfun(@(x) contains(x, 'ILMF'), {tifFiles.name}));

    for j = 1:length(tifFiles)
        imagenum = imagenum + 1;
        filePath = fullfile(subdirPath, tifFiles(j).name);

        % 读取GeoTIFF文件
        [data, R] = readgeoraster(filePath);

        % 设置0值为NaN
        data(data == 0) = NaN;

        % 更新总的地理范围
        minLat = min(minLat, R.LatitudeLimits(1));
        maxLat = max(maxLat, R.LatitudeLimits(2));
        minLon = min(minLon, R.LongitudeLimits(1));
        maxLon = max(maxLon, R.LongitudeLimits(2));

        % 将数据加入列表
        allRasters{end+1} = struct('Raster', data, 'Ref', R);
    end
end

fprintf("使用影像数量：%d\n", imagenum);


% 计算新栅格的大小
refRaster = allRasters{1}.Ref; % 使用第一个栅格作为参考
cellSizeLat = refRaster.CellExtentInLatitude;
cellSizeLon = refRaster.CellExtentInLongitude;

% 确保栅格尺寸是整数
latRange = ceil((maxLat - minLat) / cellSizeLat);
lonRange = ceil((maxLon - minLon) / cellSizeLon);

% 创建一个大的栅格容器
bigRaster = NaN(latRange, lonRange);

% 创建新的栅格引用
newRefRaster = georasterref('RasterSize', size(bigRaster), ...
    'LongitudeLimits', [minLon, maxLon], ...
    'LatitudeLimits', [minLat, maxLat]);

% 镶嵌过程
for k = 1:length(allRasters)
    raster = allRasters{k};
    [cols, rows] = geographicToIntrinsic(newRefRaster, raster.Ref.LatitudeLimits, raster.Ref.LongitudeLimits);

    % 确保索引在有效范围内
    rows = latRange - round(rows);
    cols = round(cols);
    % 稍微调整索引范围来确保一致性
    row_start = max(1, rows(2));
    row_end = min(size(bigRaster, 1), rows(1));
    col_start = max(1, cols(1));
    col_end = min(size(bigRaster, 2), cols(end));
    % 对于重叠区域，采用忽略NaN的平均
    overlap = bigRaster(row_start:row_end, col_start:col_end);
    % 确保数据维度一致
    oldvalues = double(raster.Raster);
    oldvalues(oldvalues == 0) = nan;
    newValuesResized = imresize(oldvalues, size(overlap));
    %% 异常值去除
    if(~strcmp(GasType, 'RGB'))
        newValuesResized(newValuesResized > cmax) = nan;
        newValuesResized(newValuesResized <= 0) = nan;
    end
    %%
    combined = cat(3, overlap, newValuesResized);
    bigRaster(row_start:row_end, col_start:col_end) = nanmean(combined, 3);
end
bigRaster(isnan(bigRaster)) = 0;
bigRaster_flip = flipud(bigRaster); % 上下翻转
finalRaster = bigRaster_flip;
%% 后处理
if(~strcmp(GasType, 'RGB'))
    cutoffFrequency = 8e2;
    % 应用高斯低通滤波
    filteredbigRaster = gaussianLowPassFilter(bigRaster_flip, cutoffFrequency);
    filteredbigRaster(filteredbigRaster < 0.01*cmax) = 0;

    % 8邻域版
    neighbour2 = 8;
    % 使用二值图像标记连通区域
    labeled_image = bwlabel(filteredbigRaster > 0, neighbour2); % 标记连通区域
    region_props = regionprops(labeled_image, 'Area'); % 获取各个区域的像素数目
    mask = ismember(labeled_image, find([region_props.Area] >= 100));
    maskedbigRaster = filteredbigRaster .* mask; % 筛选结果
    finalRaster = maskedbigRaster;
end
%%
figure('Position', [400 60 1200 800])
imagesc(finalRaster)
axis square
colorbar
if(strcmp(GasType, 'CO2'))
    clim([0 cmax])
elseif(strcmp(GasType, 'CH4'))
    clim([0 cmax])
end

%%
% 输出结果
outputFile = fullfile('C:\Users\huayi\Desktop\东营CCUS', [GasType, '_MergedResult.tif']);
% 现在，使用 geotiffwrite 将数据保存为 GeoTIFF
geotiffwrite(outputFile, double(finalRaster), newRefRaster);
disp('影像镶嵌完成！');
toc;