煤矸双能X射线透射图像分割方法

发表于 分类于 ， 阅读次数： 本文字数： 4.4k 阅读时长 ≈ 4 分钟

​ 煤矸双能X射线透射图像经过去噪校正后一般背景比较单一均匀，采用阈值法即可有效分割出前景目标。待分割图像如下图所示：

Snipaste_2024-04-08_20-11-34.png

​ 上图成像特点体现为明暗不均，类似于可见光成像的光照不均匀。煤密度小成像亮，贴近背景；矸石密度大成像暗，较好分割。煤矸厚度越大成像也越暗，因此阈值的选择尤为关键。既要保障分割全部目标，又要保障分割出全部的目标区域。

分割方法一：基于空皮带图像背景区域最小值的目标分割

​ 基于空皮带图像背景区域最小值的目标分割，代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
clc 
clear all;
close all;
%读取高低能图像
I = imread('C:\Users\雷子\Desktop\煤矸双能X射线透射图像分割方法\2020_8_20_15_52_11_144_H.bmp');  
I=I(:,:,1); %空皮带背景区域
I_k= I(800:1000,500:600);

figure(1),subplot(1,2,1),imshow(I);
subplot(1,2,2),imshow(I_k);

%采用中值滤波去噪
I= medfilt2(I,[3,3]); 

%计算分割阈值
min_I_k = min(min(I_k)); % min_I_k=188

%基于最小空皮带灰度值的图像分割
[m,n] = size(I);
img_I = I;
for i = 1 : m
    for j = 1 : n
        if img_I (i,j) > min_I_k 
            img_I (i,j) = 0;
        else 
             img_I (i,j) = 1;
        end
    end
end

%删除小面积区域
img_I_d= bwareaopen(img_I,10);
figure(2),subplot(1,2,1),imshow(img_I,[0,1]);
subplot(1,2,2);imshow(img_I_d);

​ 分割后的图像如下图所示：

Snipaste_2024-04-08_20-32-09.png

​ 从分割后的图像中可以看到，采用基于空皮带图像背景区域最小值的目标分割方法可以分割出全部的目标，但是目标存在有粘连区域。

分割方法二：OTSU阈值分割（大津法）

​ 采用大津法分割图像，又称OTSU阈值分割法，通过MATLAB中graythresh获取阈值，然后基于阈值thersh_H分割图像。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
clc 
clear all;
close all;
%读取高低能图像
I = imread('C:\Users\雷子\Desktop\煤矸双能X射线透射图像分割方法\2020_8_20_15_52_11_144_H.bmp');  

I=I(:,:,1);
%计算二值化图像
thersh_H = graythresh(I); %最大化类间方差 OTSU阈值分割 大津法
High_im_bw = im2bw(I,thersh_H);%thersh_H=0.4510，归一化后的阈值
%二值化取反
High_im_bw = ~High_im_bw;
High_im_bw = bwareaopen(High_im_bw,50);
imshow(High_im_bw);

%删除小面积区域
img_I_d= bwareaopen(High_im_bw,10);
imshow(img_I_d);

​ 分割后的图像如下图所示：

Snipaste_2024-04-08_20-37-01.png

​ 从分割后的图中可以看到，部分小目标由于成像接近背景而未被有效分割。

分割方法三：迭代阈值分割

​ 迭代阈值分割通过迭代方式计算最佳阈值，具有一定的自适应性，迭代阈值分割主要步骤如下：

​ ①设定初始参数T0，初始阈值T1；

​ ②使用阈值T1分割图像，G1部分像素值低于T1，G2部分像素值高于T1；

​ ③计算G1和G2各自像素平均值m1和m2，更新阈值T2=(m1+m2)/2；

​ ④若|T1-T2|<T0，则T2为最优阈值，否则令T1=T2，重复②~④。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
clc 
clear all;
close all;
%读取高低能图像
I = imread('C:\Users\雷子\Desktop\煤矸双能X射线透射图像分割方法\2020_8_20_15_52_11_144_H.bmp');  
I=im2double(I);
% 第一步:设置参数T0,并选择一个初始的估计阈值T1(取图像I像素值的最小值和最大值的平均值)
T0=0.01;
T1=(min(I(:))+max(I(:)))/2;
% 第二步：用阈值T1分割图像.将图像分成两部分:r1由灰度值大于T1的像素组成,r2是由灰度小于或等于T1的像素组成
r1=find(I>T1);% find函数返回素有非零元素的位置
r2=find(I<=T1);
% 第三步:计算r1和r2中所有像素的平均灰度值h1和h2以及新的阈值T2=(h1+h2)/2
T2=(mean(I(r1))+mean(I(r2)))/2;
% 第四步:若|T2-T1|<T0,则推出T2即为最优阈值;否则,将T2赋值给T1,并重复步骤2-4直到获取最优阈值
if abs(T2-T1)<T0
    J=imbinarize(I,T2);  % 使用imbinarize函数进行图像分割
else
    while abs(T2-T1)>=T0
        T1=T2;
        r1=find(I>T1);
        r2=find(I<=T1);
        T2=(mean(I(r1))+mean(I(r2)))/2;
    end
    J=imbinarize(I,T2);  % 使用imbinarize函数进行图像分割
end
J=~J;
subplot(121),imshow(I);
title('原始图像');
subplot(122),imshow(J);
title('迭代式阈值分割后的图像');

​ 分割结果如下如：

Snipaste_2024-04-08_20-48-14.png

​ 采用迭代阈值分割结果与采用OTSU阈值分割结果相当。

分割方法四：分块OTSU阈值分割

​ 分块OTSU阈值分割通过将图像分割成若干块，对每一块计算类间最大方差分别选取阈值来削弱成像明暗不均匀造成的漏分割。分块后需要判断分块中只有背景的情况，排除对这种分块的处理。当图像中只有背景或只有煤矸时，灰度值比较接近，采用OTSU算出的“背景”和“前景”平均灰度差Δm会很小。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
clc 
clear all;
close all;
%读取高低能图像
I = imread('C:\Users\雷子\Desktop\煤矸双能X射线透射图像分割方法\2020_8_20_15_52_11_144_H.bmp');  
I=I(:,:,1);

rn = 4;cn = 4; %定义分块区域 4X4即分成16块
[R , C] = size(I);%分别返回行和列数
rblk = R/rn;cblk = C/cn;%小块的行数和列数
x = 0:cblk:C+1;
y = 0:rblk:R+1;
M = meshgrid(x,y); %产生网格
N = meshgrid(y,x);  %产生网格
figure,imshow(I),hold on;
plot(x,N,'g'); %画出水平横线
plot(M,y,'g'); %画出垂直竖线

T = zeros(rn,cn);
dif = zeros(rn,cn);
J = false(R,C);%初始化二值图
X = uint8(zeros(rblk,cblk));

%分块阈值，并判断类间灰度差以排除纯背景或纯物体的干扰
for r=1:rn
    for c=1:cn
        r0=rblk*(r-1)+1;r1=rblk*r;
        c0=cblk*(c-1)+1;c1=cblk*c;
        X = I(r0:r1,c0:c1);
        T(r,c) = graythresh(X);                             %计算分割阈值
        [h,~] = histcounts(X,0:255);
        T_int =uint8(T(r,c)*255);
        dif(r,c) = graydiffer(h,T_int);                     %计算类间灰度差的函数略
        str = ['T=',num2str(T_int),' △m=',num2str(dif(r,c))];
        text(c0+cblk/4,r0+rblk/2,str,'color','black');%显示信息
        if dif(r,c)>30 %这里设定m
            J(r0:r1,c0:c1) = ~im2bw(X,T(r,c));
        end
    end
end

bw_F = bwareaopen(J,50); %删除小面积连通域区域
bw_F  = J;
figure,imshow(bw_F); %显示最终分割结果

​ 注意！调用函数graythresh：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
%功能：计算一幅图像前景和背景类间平均灰度差
%输入：直方图数据h,分割阈值T
%输出：类间平均灰度差
%作者：LEI HE
%日期：2024/04/08
function[differ] = graydiffer(h,T)
    s1 = sum(h(1:T));
    s2 = sum(h(T:255));
    n1 = 1:T;
    n2 = T:255;
    u1 = double(n1)*h(1:T)' / s1; %背景灰度均值
    u2 = double(n2)*h(T:255)' / s2; %前景灰度均值
    differ = uint8(u2-u1);
end

​ 分块图像如下图所示：

Snipaste_2024-04-08_20-59-00.png

​ 分块OTSU阈值分割算法将图像分为四列，保证每一个块中都包含背景区域，在此基础上降低分块数量可提高计算效率。图中T为单个区域分块的最佳分割阈值，Δm为此分块对应类间灰度差。设定类间灰度差Δm为一个合适的阈值，如Δm=30，即只对Δm>30的分块进行阈值二值化处理，而Δm<30的分块全部像素点均置零。

​ 分割后的图像如下图所示：

Snipaste_2024-04-08_20-55-02.png

​ 采用分块OTSU阈值分割，在保证分割全部目标的同时，不会引入大量的粘连。

​ 原图像的灰度直方图：

Snipaste_2024-04-08_21-13-35.png

​ 使用背景阈值分割将分布在150附近的灰度值划分到前景中，这是粘连的主要原因。采用背景阈值分割获得的区域边缘一定为目标X射线成像的真实边缘，而分块OTSU阈值分割会舍弃此部分边缘。X射线透射成像在不遇到物体时背景值不会降低，一但像素值低于背景像素值那么则说明有物体被穿透。但为了避免引入大量不必要的重叠粘连，选择分块OTSU阈值分割对图像进行二值化处理。

1.png

​ 本文源码对应博士论文第三章，第5节内容。