function [TIFF, flag] = MaskPlume2(TIFF)
DISTANCE = 0;
neighbour = 8; % 使用8邻域处理
distances = bwdist(TIFF.maskdata > 0); % 计算到每个非零像素的距离
connected = distances <= DISTANCE; % 距离小于某一数值的像素被认为是连通的
% 使用二值图像标记连通区域
labeled_image = bwlabel(connected, neighbour); % 标记连通区域
%% 提取主羽流
% 默认羽流位置为(TIFF.rowindex, TIFF.colindex)
if isfield(TIFF, 'rowIndex_new') && isfield(TIFF, 'colIndex_new')
    x0 = TIFF.colIndex_new;
    y0 = TIFF.rowIndex_new;
else
    x0 = TIFF.colIndex;
    y0 = TIFF.rowIndex;
end
% 计算每个点到指定点的欧氏距离
[rows, cols] = size(TIFF.maskdata);
[xx, yy] = meshgrid(1:cols, 1:rows); % 生成与图像大小相同的网格
distances_to_facility = sqrt((yy - y0).^2 + (xx - x0).^2); % 计算每个点到设施点的距离

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)) > 30)
        continue; % 如果该区域到设施的最近距离超过30格(900m)，跳过
    end
    temp_maskdata = TIFF.maskdata(region_mask); % 当前区域的数据
    region_weighted_strength = sum(temp_maskdata(:)) / mean(distances_to_facility(region_mask) + 1)^1.5;
    if (region_weighted_strength > max_strength) % 如果超过当前最大强度
        max_strength = region_weighted_strength; % 记录强度
        main_plume_mask = region_mask; % 记录掩膜
        main_plume_label = label; % 记录编号
    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(:));
[y_max, x_max] = ind2sub(size(main_plume_mask), max_index);
%% 主羽流信息
if(exist('main_plume_label', 'var'))
    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);
    %% 以投影点为起点，投影点到羽流中心的向量为方向，定义射线
    if(norm([x_max, y_max] - [x0, y0]) < 50) % 如果纠正的设施点与原来相差不大
        segment_start = [x_max, y_max];
    else % 否则采用设施点
        segment_start = [x0, y0];
    end
    parameter = 3; % 控制线段比
    segment_end =  main_plume_centroid + parameter * (main_plume_centroid - segment_start);
    vector_main_plume = segment_end - segment_start;
    segment_length = norm(vector_main_plume);
    %% 提取次羽流
    % 初始值
    sub_plume_mask = main_plume_mask;
    sub_plume_label = main_plume_label;
    for label = 1:max(labeled_image(:))
        region_mask = labeled_image == label;
        % 计算当前区域到主羽流的最短距离
        distance_map = bwdist(main_plume_mask) .* region_mask;
        distance_map(distance_map == 0) = nan;
        distance_to_main_plume = min(distance_map(:)); % 找到最短距离
        if(distance_to_main_plume > 30) % 当前羽流到主羽流的最短距离不得超过30
            continue;
        end
        if(label == main_plume_label) % 如果是主羽流
            continue;
        end
        plume_distance_to_segment_min = region_to_segment_distance(region_mask, segment_start, segment_end);
        % 如果当前区域到主羽流线段的距离小于线段一半长度
        % 并且区域中心与排放设施的构成的向量，排放设施与主羽流的构成的向量，两者夹角小于60度即可
        current_centroid = region_centroid(label).Centroid;
        vector_current_plume = current_centroid - segment_start;
        angle = vector_angle(vector_main_plume, vector_current_plume);
        
        if(segment_length < 80)
            thresh_angle = pi/6;
        elseif(segment_length < 160)
            thresh_angle = pi/8;
        else
            thresh_angle = pi/12;
        end
        if( plume_distance_to_segment_min < parameter * segment_length && angle < thresh_angle )
            sub_plume_mask = region_mask | sub_plume_mask;
        end
    end
    % 将最大和次大的羽流合并
    combined_mask = main_plume_mask | sub_plume_mask;
    % 使用掩模来访问原始图像的区域
    sub_plume_centroid = region_centroid(sub_plume_label).Centroid;
    x2 = sub_plume_centroid(1);
    y2 = sub_plume_centroid(2);
    %% 膨胀区域
    % 定义膨胀的结构元素（例如，一个矩形）
    se = strel('rectangle', [3, 3]); % 3x3 的矩形
    % 对二值图像执行膨胀操作
    dilated_mask = imdilate(combined_mask, se);
    
    %% 填充孔洞
    filled_mask = imfill(dilated_mask, 'holes');
    filled_plume_data = TIFF.maskdata;
    filled_plume_data(~filled_mask) = 0; % 将不在最近连通区域内的像素设为0
    
    %% 赋值
    TIFF.finalmask = filled_plume_data > 0;
    TIFF.finaldatashow = filled_plume_data;
    TIFF.finaldata = TIFF.cutdata .* TIFF.finalmask;
    TIFF.x1 = x1;
    TIFF.y1 = y1;
    TIFF.x2 = x2;
    TIFF.y2 = y2;
    TIFF.xmax = x_max;
    TIFF.ymax = y_max;
    TIFF.segment_start = segment_start;
    TIFF.segment_end = segment_end;
    %%
    flag = true;
    %% 计算反演不确定度
    bkgdata = TIFF.cutdata(~TIFF.finalmask); % 剔除羽流部分
    bkgdata(bkgdata == 0) = [];
    uncertainty = std(bkgdata(:));
    TIFF.finaldata_u = uncertainty / 0.6120; % 因为截取了负半轴

else
    flag = false;
    TIFF = -1;
end
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
%%
function min_distance = region_to_segment_distance(region, segment_start, segment_end)
% 初始化最小距离为正无穷
min_distance = Inf;

% 遍历区域内的每个像素点
[rows, cols] = size(region);
for row = 1:rows
    for col = 1:cols
        % 如果该像素属于区域
        if region(row, col)
            % 计算当前像素点到线段的距离
            point = [col, row];
            distance = distance_to_segment(point, segment_start, segment_end);
            
            % 更新最小距离
            min_distance = min(min_distance, distance);
        end
    end
end
end

%% 计算向量夹角
function angle = vector_angle(vector1, vector2)
% 计算向量的点积
dot_product = dot(vector1, vector2);

% 计算向量的模长
norm_vector1 = norm(vector1);
norm_vector2 = norm(vector2);

% 计算夹角的余弦值
cos_angle = dot_product / (norm_vector1 * norm_vector2);

% 使用反余弦函数计算夹角的弧度值
angle = acos(cos_angle);
end