function TIFF = MaskPlume(TIFF, thresh)
%% 二值提取
% DISTANCE 距离小于该数值的像素被认为是连通的
DISTANCE = 1;
neighbour = 4; % 使用四邻域处理
distances = bwdist(TIFF.threshdata > 0); % 计算到每个非零像素的距离
connected = distances <= DISTANCE; % 距离小于某一数值的像素被认为是连通的
% 使用二值图像标记连通区域
labeled_image = bwlabel(connected, neighbour); % 使用 4-邻域连接法则标记连通区域
region_props = regionprops(labeled_image, 'Area'); % 获取各个区域的像素数目
% 仅保留像素数目达到 thresh 个的部分
TIFF.mask = ismember(labeled_image, find([region_props.Area] >= thresh));
while(sum(TIFF.mask(:)) == 0) % 如果没有符合条件的
    thresh = thresh - 4;
    TIFF.mask = ismember(labeled_image, find([region_props.Area] >= thresh));
    if(thresh <= 0)
        break;
    end
end
TIFF.maskdata = TIFF.threshdata .* TIFF.mask; % 筛选结果
%% 提取主羽流
% 假设指定位置为(TIFF.rowindex, TIFF.colindex)
x0 = TIFF.rowIndex;
y0 = TIFF.colIndex;
% 计算每个点到指定点的欧氏距离
[rows, cols] = size(TIFF.cutdata);
[xx, yy] = meshgrid(1:cols, 1:rows); % 生成与图像大小相同的网格
distances_to_facility = sqrt((yy - x0).^2 + (xx - y0).^2); % 计算每个点到设施点的距离
% 获取设施点所属的连通区域标签
facility_label = labeled_image(x0, y0);
% 如果指定点属于某个连通区域
if facility_label ~= 0
    main_plume_mask = labeled_image == facility_label;
    main_plume_label = facility_label;
else % 找到最近最强的连通区域
    max_strength = 0; % 初始值
    main_plume_mask = false(size(labeled_image));
    for label = 1 : max(labeled_image(:))
        region_mask = labeled_image == label; % 当前区域的掩膜
        % 找到当前区域中像素的最大强度及其坐标
        if(min(distances_to_facility(region_mask)) > 120)
            continue; % 如果该区域到设施的最近距离超过120，跳过
        end
        temp_maskdata = TIFF.maskdata(region_mask); % 当前区域的数据
        region_weighted_strength = sum(temp_maskdata(:)) / mean(distances_to_facility(region_mask));
        if (region_weighted_strength > max_strength) % 如果超过当前最大强度
            max_strength = region_weighted_strength; % 记录强度
            main_plume_mask = region_mask; % 记录掩膜
            main_plume_label = label; % 记录编号
        end
    end
end
% 获取主羽流排放数值
main_plume_data = TIFF.maskdata;
main_plume_data(~main_plume_mask) = 0; % 将不在最近连通区域内的像素设为0
% 获取主羽流到排放设施的距离
main_plume_distance = distances_to_facility;
main_plume_distance(~main_plume_mask) = 9999; % 非羽流设置为高值
%% 根据主羽流修正排放源位置
% 找到当前区域中像素的最大最近强度及其坐标->用于修正排放设施位置
[~, max_index] = max(main_plume_data(:) - main_plume_distance(:));
[y_max, x_max] = ind2sub(size(main_plume_mask), max_index);
%% 主羽流信息
region_centroid = regionprops(labeled_image, 'Centroid'); % 获取羽流的中心信息
main_plume_centroid = region_centroid(main_plume_label).Centroid; % 主羽流的中心
x1 = main_plume_centroid(1);
y1 = main_plume_centroid(2);
% region_orientation = regionprops(labeled_image, 'Orientation');
% first_orientation = region_orientation(nearest_label).Orientation;
% theta = deg2rad(first_orientation); % 将角度转换为弧度
% 计算坐标位置到第一羽流方向直线的投影点
% L = @(x, y) x*sin(theta) - y*cos(theta) - (x1*sin(theta) - y1*cos(theta));
% 计算投影点坐标
% x_projection = (y0 - y1)*sin(theta)*cos(theta)+x0*cos(theta)^2+x1*sin(theta)^2;
% y_projection = y0*sin(theta)^2 + y1*cos(theta)^2 + (x0-x1)*sin(theta)*cos(theta);
%% 以投影点为起点，投影点到羽流中心的向量为方向，定义射线
ray_start = [x_max, y_max];
ray_direction = main_plume_centroid - ray_start;

%% 提取次羽流
% 初始值
max_strength_ratio = 0;
sub_plume_mask = main_plume_mask;
sub_plume_label = main_plume_label;
if(sum(main_plume_mask(:)) < 2000) % 如果主羽流不够强，需要提取次羽流
    for label = 1:max(labeled_image(:))
        region_mask = labeled_image == label;
        % 排除主羽流块
        if isequal(region_mask, main_plume_mask)
            continue;
        end
        % 计算当前羽流块的中心到直线L的距离
        current_centroid = region_centroid(label).Centroid;
%         distance_to_main_plume = norm(current_centroid - main_plume_centroid); % 计算到第一羽流中心的距离
        
        % 计算当前区域到主羽流的最短距离
        distance_map = bwdist(main_plume_mask) .* region_mask;
        distance_to_main_plume = min(distance_map(:)); % 找到最短距离
        if(distance_to_main_plume > 20) % 当前羽流到主羽流的最短距离不得超过20
            continue;
        end
        % 计算羽流块的总强度
        current_plume_data = TIFF.maskdata(region_mask);
        region_strength = sum(current_plume_data(:));
        plume_distance_to_ray = distance_to_ray(current_centroid, ray_start, ray_direction); % 调用函数计算点到射线的距离

        % 计算强度与距离比
        strength_ratio = region_strength^2 / (plume_distance_to_ray + distance_to_main_plume);
        
        % 更新最大比率和最强的羽流块
        if strength_ratio > max_strength_ratio
            max_strength_ratio = strength_ratio;
            sub_plume_mask = region_mask;
            sub_plume_label = label;
        end
    end
    % 将最大和次大的羽流合并
    combined_mask = main_plume_mask | sub_plume_mask;
    % 使用掩模来访问原始图像的区域
    combined_plume_data = TIFF.maskdata;
    combined_plume_data(~combined_mask) = 0; % 将不在最近连通区域内的像素设为0
else
    combined_mask = main_plume_mask; % 整合羽流掩膜
    combined_plume_data = main_plume_data; % 整合羽流数据
end
sub_plume_centroid = region_centroid(sub_plume_label).Centroid;
x2 = sub_plume_centroid(1);
y2 = sub_plume_centroid(2);
%% 赋值
TIFF.finaldata = combined_plume_data;
TIFF.finalmask = combined_mask;
TIFF.x1 = x1;
TIFF.y1 = y1;
TIFF.x2 = x2;
TIFF.y2 = y2;
TIFF.xmax = x_max;
TIFF.ymax = y_max;
% TIFF.theta = theta;
TIFF.ray_start = ray_start;
TIFF.ray_direction = ray_direction;
% TIFF.x_projection = x_projection;
% TIFF.y_projection = y_projection;
%% 不确定度
uncertainty = std(TIFF.finaldata(:));
finaldata_u = ones(size(TIFF.finaldata)) * uncertainty;
TIFF.finaldata_u = finaldata_u;
end

%% 计算点到射线的距离
function distance = distance_to_ray(point, ray_start, ray_direction)
% 计算射线方向的单位向量
ray_direction_unit = ray_direction / norm(ray_direction);

% 计算投影点到射线上的参数 t
t = dot(point - ray_start, ray_direction_unit);

if t >= 0
    % 计算投影点
    projected_point = ray_start + t * ray_direction_unit;
else
    % 当投影点在射线的起点之前时，返回点到起点的距离
    projected_point = ray_start;
end

% 计算点到射线的距离
distance = norm(point - projected_point);
end

%% 计算点到线段的距离
function distance = distance_to_segment(point, segment_start, segment_end)
    % 计算线段的方向向量和长度
    segment_direction = segment_end - segment_start;
    segment_length = norm(segment_direction);
    segment_direction_unit = segment_direction / segment_length;
    
    % 计算投影点到线段起点的向量
    point_to_start = point - segment_start;
    
    % 计算投影点在线段方向上的投影长度
    t = dot(point_to_start, segment_direction_unit);
    
    if t <= 0
        % 如果投影点在线段起点之前，则最近距离是点到线段起点的距离
        distance = norm(point - segment_start);
    elseif t >= segment_length
        % 如果投影点在线段终点之后，则最近距离是点到线段终点的距离
        distance = norm(point - segment_end);
    else
        % 如果投影点在线段上，则最近距离是点到投影点的距离
        projected_point = segment_start + t * segment_direction_unit;
        distance = norm(point - projected_point);
    end
end