提取R值均值、灰度均值及对R值图像进行分割计算像素面积占比

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

本节主要内容为提取高低能图像目标区域R值均值、灰度均值等信息，并绘制R值分布图：

Snipaste_2024-05-08_10-09-31.png

首先是将R值均值、灰度均值进行绘图的一个代码：

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
45
46
47
48
49
50
51
52
53
54
clc;
clear all;
close all;
load('D:\r_hl_hh');
figure,
plot(r_hl_hh(1:35,1),'r'),hold on,
plot([35:1:70],r_hl_hh(35:70,1),'g'),hold on,
plot([70:1:100],r_hl_hh(70:100,1),'b'),hold on,

plot(r_hl_hh(:,2),'k'),hold on,
legend('Rich coal','Charred coal','Gas coal','Gangue','Thrashold-1','Thrashold-2');

% plot([0,100],[1.222215 1.222215],'k-');
plot([0,100],[1.42054 1.42054],'k-');
% plot([0,100],[1.618934 1.618934],'k-');
plot([0,100],[1.36027 1.36027],'k-');   
plot([0,100],[1.38815 1.38815],'g-');%中位数均值
plot([0,100],[1.3 1.3],'g-');
% median(r_hl_hh(:,2))
min(r_hl_hh(:,2))
% scatter(r_hl_hh(:,3),r_hl_hh(:,1)),hold on;
% scatter(r_hl_hh(:,4),r_hl_hh(:,2)),hold on;
legend('coal r value','gangue r value');

%低能灰度均值分布,66.4380
figure,
plot(r_hl_hh(1:35,3),'r'),hold on,
plot([35:1:70],r_hl_hh(35:70,3),'g'),hold on,
plot([70:1:100],r_hl_hh(70:100,3),'b'),hold on,
plot(r_hl_hh(:,4),'k'),hold on,

mean_c = mean(r_hl_hh(:,3));
mean_g = mean(r_hl_hh(:,4));
value_h = (mean_c +mean_g)/2;
plot([0,100],[value_h value_h],'k-');  
legend('Rich coal','Charred coal','Gas coal','Gangue','Thrashold');
%plot([0,100],[(min(r_hl_hh(:,3))+max(r_hl_hh(:,4)))/2 (min(r_hl_hh(:,3))+max(r_hl_hh(:,4)))/2],'k-');  
% (median(r_hl_hh(:,3)) + median(r_hl_hh(:,4)))/2;
% plot([0,100],[66.5216 66.5216],'k-');  


%高能灰度均值分布48.7378
figure,
plot(r_hl_hh(1:35,5),'r'),hold on,
plot([35:1:70],r_hl_hh(35:70,5),'g'),hold on,
plot([70:1:100],r_hl_hh(70:100,5),'b'),hold on,
plot(r_hl_hh(:,6),'k'),hold on,

mean_cc = mean(r_hl_hh(:,5));
mean_gg = mean(r_hl_hh(:,6));
value_hh = (mean_cc +mean_gg)/2;
plot([0,100],[value_hh value_hh],'k-');  
legend('Rich coal','Charred coal','Gas coal','Gangue','Thrashold');

输出图像：

QQ20240506205024.png

上面图像的绘图特点在于将多组数据绘制在一个figure里面。

接下来是R只均值和灰度均值提取代码：R.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
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
clc 
clear all;
%读取高低能图像
High_im = imread('F:\0820测试\0-30g\2020_8_19_14_22_19_915_H.bmp');
Low_im = imread('F:\0820测试\0-30g\2020_8_19_14_22_19_899_L.bmp');

%去除探测器损坏区域图像
High_im = High_im(40:1000,40:1000);
Low_im = Low_im(40:1000,40:1000);
% High_im = High_im(1:500,100:600);
% Low_im = Low_im(1:500,100:600);
figure(1),subplot(1,2,1),imshow(High_im);
subplot(1,2,2),imshow(Low_im);
%中值滤波
% High_im = medfilt2(High_im,[3,3]);
% Low_im = medfilt2(Low_im,[3,3]);
%计算二值化图像
thersh_H = graythresh(High_im);
thersh_L = graythresh(Low_im);
High_im_bw = im2bw(High_im,thersh_H);
Low_im_bw = im2bw(Low_im,thersh_L);

%二值化取反
High_im_bw = ~High_im_bw;
Low_im_bw = ~Low_im_bw ;
High_im_bw = bwareaopen(High_im_bw,50);
Low_im_bw = bwareaopen(Low_im_bw,50);
figure(2),subplot(1,2,1),imshow(High_im_bw);
subplot(1,2,2),imshow(Low_im_bw);

Low_im_nobord = imclearborder(Low_im_bw,4);
seD = strel('diamond',1);
BWfinal = imerode(Low_im_nobord,seD);
figure(3),imshow(BWfinal);
%空皮带平均灰度值构建矩阵
% em_H = double(High_im.*(0)) +77.3099; 
% em_L = double(Low_im.*(0)) +122.7394;
% 
em_H = double(High_im.*(0)) +194.4723; 
em_L = double(Low_im.*(0)) +193.9575;
% em_H = double(High_im.*(0)) +194.4812;
% em_L = double(Low_im.*(0)) +194.0422;
%计算并显示R值图像
R_im=(log(double(em_L) ./ double(Low_im)) ) ./ (log(double(em_H) ./ double(High_im)));
figure(4),imshow(R_im);
%掩膜处理提取目标区域R值图像
R_im_f = R_im .* BWfinal;
figure(5),imshow(R_im_f);

R_im_f(find(isnan(R_im_f)==1)) = 0;
R_im_f(find(isinf(R_im_f)==1)) = 0;

% R_im_f = uint8(R_im_f);
% figure,imhist(R_im_f);
figure,imagesc(uint8(R_im_f),[0,2]);colorbar
colormap jet ;

% %映射R值图像到0~255灰度区间内
% Omax = 255;
% Omin = 0;
% x = find(R_im_f <1 & R_im_f >0.8 );
% y = find(R_im_f <2 & R_im_f >0);
% p = R_im_f (x);
% q = R_im_f (y);
% a =min(p);
% b =max(q);
% for i = 40 : 1000
%     for j = 40 : 1000
%         R_im_new(i,j) =(Omax - Omin) * (R_im_f(i,j) - a) / (b - a) + Omin;
%         if R_im_new(i,j) == 0
%             R_im_new(i,j) =255;
%         end
% %         if R_im_new(i,j) > 255
% %             R_im_new(i,j) =0;
% %         end
%     end
% end      
% figure(6),imshow(uint8(R_im_new),[0,255]);
% 
% figure,imagesc(uint8(R_im_new),[0,255]);colorbar
% colormap jet ;

% R_im_new = uint8(R_im_new);
% figure,imhist(R_im_new);
%%

% %计算目标区域R值均值
[m,n] = size(R_im_f);
% I1 = R_im_f;
% for i = 1:m
%     for j = 1:n
%         if I1(i,j) > 2
%             I1(i,j) = 0;
%         end
%     end
% end
% figure(12),imshow(I1);
%计算I1R值均值,首先获取单个目标区域，背景为0，目标区域为1
I2 = bwlabel(Low_im_bw);
r = max(max(I2));
S = cell(r,1);
for w = 1:r
    I2 = bwlabel(Low_im_bw);
    for i = 1:m
        for j =1: n
            if I2(i,j) ~= w
                I2(i,j) = 0;
            end
            if I2(i,j) ==w
                I2(i,j) = 1;
            end
        end
    end
    S{w,1} = I2;
end

mean_r = cell(r,1);
for w = 1:r
    mean_r{w,1} = R_im_f .* S{w,1};
end
% figure,imshow(S{2,1});
% figure,imshow( mean_r{2,1} );

for w = 1:r
    mean_i(w,1)= sum(sum(mean_r{w,1})) / nnz(mean_r{w,1});
end

%灰度均值计算

%低能
for w = 1:r
    mean_hl{w,1} =double(Low_im) .*S{w,1};
end

for w = 1:r
    mean_hll(w,1)= sum(sum(mean_hl{w,1}))/nnz(mean_hl{w,1});
end    

%高能
for w = 1:r
    mean_hh{w,1} =double(High_im) .*S{w,1};
end

for w = 1:r
    mean_hhh(w,1)= sum(sum(mean_hh{w,1}))/nnz(mean_hh{w,1});
end 

%%
% scatter(mean_i(16:37,1),mean_hh(16:37,1)),hold on;
% scatter(mean_i(43:78,1),mean_hh(43:78,1));
% 
% %绘制三维图
% u = 1:1:20;
% v = 1:1:20;
% [x,y] = meshgrid(u,v);
% [~,z]= meshgrid(mean_i(1:20,1),mean_i(1:20,1));
% surf(x,y,z),hold on;
% [~,t] = meshgrid(mean_i(41:60,1),mean_i(41:60,1));
% surf(x,y,t);
%%
%统计各区域红色像素占比
%二值化取反图像连通域检测
[G,num] = bwlabel(Low_im_bw);
%连通域分离
EZ = cell(num,1);
for c = 1:num
    G = bwlabel(Low_im_bw,8);
    for i = 1:m
        for j = 1:n
             if G(i,j) ~= c
                 G(i,j) = 0;  
             end
        end
    end
    EZ{c,1} = G .* (1/c);
end

%imshow( EZ{1,1})
HL = cell(num,1);
for c = 1:num
    HL{c,1} = J9.* EZ{c,1};
end

%figure(3),imshow(HL{18,1},[])
%计算红色像素占比
HL_L = cell(num,1);
for c = 1:num
    p = 0;
    q = 0;
    for i = 1:m
        for j = 1:n
            if  HL{c,1} (i,j) == 1
              p = p +1  ;
            end
            if HL{c,1}(i,j) == 2
                q = q+1;
            end
        end
    end
    HL_L{c,1} = q/(p+q);
end
for c= 1:num
    GS(c,1) = HL_L{c,1};
end

接下来是通过设定阈值来分割R值图像，然后计算不同分割阈值下，单个目标中不同阈值区间像素点占比，R_fg.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
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
clc 
clear all;
%读取高低能图像
High_im = imread('D:\2019_1_16_16_27_42_970_H.bmp');
Low_im = imread('D:\2019_1_16_16_27_42_970_L.bmp');
High_im = High_im(1:1000,1:550);
Low_im = Low_im(1:1000,1:550);
figure,
subplot(1,2,1),imshow(High_im);
subplot(1,2,2),imshow(Low_im);
%中值滤波
High_im = medfilt2(High_im,[3,3]);
Low_im = medfilt2(Low_im,[3,3]);
%计算二值化图像
thersh_H = graythresh(High_im);
thersh_L = graythresh(Low_im);
High_im_bw = im2bw(High_im,thersh_H);
Low_im_bw = im2bw(Low_im,thersh_L);
%二值化取反
High_im_bw = ~High_im_bw;
Low_im_bw = ~Low_im_bw ;
% High_im_bw = bwareaopen(High_im_bw,5);
% Low_im_bw = bwareaopen(Low_im_bw,5);
figure,
subplot(1,2,1),imshow(High_im_bw);
subplot(1,2,2),imshow(Low_im_bw);

% Low_im_nobord = imclearborder(Low_im_bw,4);
% seD = strel('diamond',1);
% BWfinal = imerode(Low_im_nobord,seD);
% figure,
%  imshow(BWfinal);
%空皮带平均灰度值构建矩阵
% em_im_H = imread('D:\2020_7_31_9_30_26_154_H.bmp');
% em_im_L= imread('D:\2020_7_31_9_30_26_154_L.bmp');
% em_im_H = em_im_H(300:300,400:400);
% mean_H = mean(mean(em_im_H));
% em_im_L = em_im_L(300:300,400:400);
% mean_L = mean(mean(em_im_L));
em_H = double(High_im.*(0)) +77.3099;
em_L = double(Low_im.*(0)) +122.7394;
%计算并显示R值图像
R_im= (log(double(em_L) ./ double(Low_im)) )./ (log(double(em_H) ./ double(High_im)));
figure,
imshow(R_im);
%掩膜处理提取目标区域R值图像
R_im_f = R_im .* Low_im_bw;
figure,
imshow(R_im_f);
% R_im_f = uint8(R_im_f);
% figure,imhist(R_im_f);
%将R值图像中NaN置零
R_im_f(find(isnan(R_im_f)==1)) = 0;

%阈值分割
img = R_im_f;
[M,N]=size(img);
out=zeros(M,N);
 
% t1 =1.210232401;   %mean_coal_R
% t2 = 1.64980168;  %mean_gangue_R
% t = (t1 + t2 )/2;
t1 = 1.3
t2 = 1.42
 for m=1:M
     for n=1:N
         if (img(m,n)<t1 )
             out(m,n)=1;
         end
         if(img(m,n)>t1 && img(m,n)<t2 )
                 out(m,n)=3;  %灰度级可调
         end
         if(img(m,n) == 0)
                 out(m,n)=0;  %灰度级可调
         end
         if(img(m,n)>t2)
                 out(m,n)=2;  %灰度级可调
         end
      end
 end
figure,imshow(out);
 %伪色彩处理
I=double(out);
[m,n]=size(I);
R = zeros(m,n);
G = zeros(m,n);
B = zeros(m,n);
for i=1:m
    for j=1:n
        %赋予背景层白色
               if I(i,j) ==0
                    R(i,j)=255;
                    G(i,j)=255;
                    B(i,j)=255;
               end
         %赋予煤层绿色
               if I(i,j)==1
                  R(i,j)=0;
                  G(i,j)=255;
                  B(i,j)=0;
               end
         %赋予矸石层红色             
               if I(i,j)==2
                  R(i,j)=255;
                  G(i,j)=0;
                  B(i,j)=0;
               end
          %赋予中间层蓝色             
               if I(i,j)==3
                  R(i,j)=0;
                  G(i,j)=0;
                  B(i,j)=255;
               end
           end      
end

rgbim = zeros(m,n,3);
for i=1:m
    for j=1:n
        rgbim(i,j,1)=R(i,j);
        rgbim(i,j,2)=G(i,j);
        rgbim(i,j,3)=B(i,j);
    end
end

rgbim=rgbim / 256;
subplot(1,2,1);imshow(High_im);
subplot(1,2,2);
imshow(rgbim);
hold on ;

%统计各区域红色像素占比
%二值化取反图像连通域检测
[G,num] = bwlabel(Low_im_bw);
%连通域分离
EZ = cell(num,1);
for c = 1:num
    G = bwlabel(Low_im_bw,8);
    for i = 1:m
        for j = 1:n
             if G(i,j) ~= c
                 G(i,j) = 0;  
             end
        end
    end
    EZ{c,1} = G .* (1/c);
end
%imshow( EZ{1,1})
HL = cell(num,1);

for c = 1:num
    HL{c,1} = out.* EZ{c,1};
end
%figure,imshow(HL{21,1},[])
%计算红色像素占比
HL_L = cell(num,1);
for c = 1:num
    p = 0;
    q = 0;
    for i = 1:m
        for j = 1:n
            if  HL{c,1} (i,j) == 1
              p = p +1  ;
            end
            if HL{c,1}(i,j) == 2
                q = q+1;
            end
        end
    end
    HL_L{c,1} = q/(p+q);
end
for c= 1:num
    GS(c,1) = HL_L{c,1};
end
%矸石含量可视化
plot(GS)

R_yxfg.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
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
clc 
clear all;
%读取高低能图像
% High_im = imread('D:\双能X射线\煤(20-50mm)\2020_7_31_9_22_13_928_H.bmp');
% Low_im = imread('D:\双能X射线\煤(20-50mm)\2020_7_31_9_22_13_928_L.bmp');
High_im = imread('D:\双能X射线\煤夹矸(20-50mm)\2020_7_31_9_30_37_947_H.bmp');
Low_im = imread('D:\双能X射线\煤夹矸(20-50mm)\2020_7_31_9_30_37_947_L.bmp');
High_im = High_im(40:1000,80:1000);
Low_im = Low_im(40:1000,80:1000);
figure,
subplot(1,2,1),imshow(High_im);
subplot(1,2,2),imshow(Low_im);
%中值滤波
High_im = medfilt2(High_im,[3,3]);
Low_im = medfilt2(Low_im,[3,3]);
%计算二值化图像
thersh_H = graythresh(High_im);
thersh_L = graythresh(Low_im);
High_im_bw = im2bw(High_im,thersh_H);
Low_im_bw = im2bw(Low_im,thersh_L);
%二值化取反
High_im_bw = ~High_im_bw;
Low_im_bw = ~Low_im_bw ;
% High_im_bw = bwareaopen(High_im_bw,5);
% Low_im_bw = bwareaopen(Low_im_bw,5);
figure,
subplot(1,2,1),imshow(High_im_bw);
subplot(1,2,2),imshow(Low_im_bw);

% Low_im_nobord = imclearborder(Low_im_bw,4);
% seD = strel('diamond',1);
% BWfinal = imerode(Low_im_nobord,seD);
% figure,
%  imshow(BWfinal);
%空皮带平均灰度值构建矩阵
% em_im_H = imread('D:\2020_7_31_9_30_26_154_H.bmp');
% em_im_L= imread('D:\2020_7_31_9_30_26_154_L.bmp');
% em_im_H = em_im_H(300:300,400:400);
% mean_H = mean(mean(em_im_H));
% em_im_L = em_im_L(300:300,400:400);
% mean_L = mean(mean(em_im_L));
em_H = double(High_im.*(0)) +191;
em_L = double(Low_im.*(0)) +195;
%计算并显示R值图像
R_im= (log(double(em_L) ./ double(Low_im)) )./ (log(double(em_H) ./ double(High_im)));
figure,
imshow(R_im);
%掩膜处理提取目标区域R值图像
R_im_f = R_im .* Low_im_bw;
figure,
imshow(R_im_f);
% R_im_f = uint8(R_im_f);
% figure,imhist(R_im_f);
%将R值图像中NaN置零
R_im_f(find(isnan(R_im_f)==1)) = 0;

%阈值分割
img = R_im_f;
[M,N]=size(img);
out=zeros(M,N);
 
t1 =1.3357;   %mean_coal_R
t2 = 1.6003;  %mean_gangue_R
t = (t1 + t2 )/2;
 for m=1:M
     for n=1:N
         if (img(m,n)<t && img(m,n) ~= 0)
             out(m,n)=1;
         end
         if(img(m,n)>t)
                 out(m,n)=2;  %灰度级可调
         end
         if(img(m,n) == 0)
                 out(m,n)=0;  %灰度级可调
         end
      end
 end
figure,imshow(out);
 %伪色彩处理
I=double(out);
[m,n]=size(I);
R = zeros(m,n);
G = zeros(m,n);
B = zeros(m,n);
for i=1:m
    for j=1:n
        %赋予背景层白色
               if I(i,j) ==0
                    R(i,j)=255;
                    G(i,j)=255;
                    B(i,j)=255;
               end
         %赋予煤层绿色
               if I(i,j)==1
                  R(i,j)=0;
                  G(i,j)=255;
                  B(i,j)=0;
               end
         %赋予矸石层红色             
               if I(i,j)==2
                  R(i,j)=255;
                  G(i,j)=0;
                  B(i,j)=0;
               end
           end      
end

rgbim = zeros(m,n,3);
for i=1:m
    for j=1:n
        rgbim(i,j,1)=R(i,j);
        rgbim(i,j,2)=G(i,j);
        rgbim(i,j,3)=B(i,j);
    end
end

rgbim=rgbim / 256;
subplot(1,2,1);imshow(High_im);
subplot(1,2,2);
imshow(rgbim);
hold on ;

%统计各区域红色像素占比
%二值化取反图像连通域检测
[G,num] = bwlabel(Low_im_bw);
%连通域分离
EZ = cell(num,1);
for c = 1:num
    G = bwlabel(Low_im_bw,8);
    for i = 1:m
        for j = 1:n
             if G(i,j) ~= c
                 G(i,j) = 0;  
             end
        end
    end
    EZ{c,1} = G .* (1/c);
end
%imshow( EZ{1,1})
HL = cell(num,1);

for c = 1:num
    HL{c,1} = out.* EZ{c,1};
end
%figure,imshow(HL{21,1},[])
%计算红色像素占比
HL_L = cell(num,1);
for c = 1:num
    p = 0;
    q = 0;
    for i = 1:m
        for j = 1:n
            if  HL{c,1} (i,j) == 1
              p = p +1  ;
            end
            if HL{c,1}(i,j) == 2
                q = q+1;
            end
        end
    end
    HL_L{c,1} = q/(p+q);
end
for c= 1:num
    GS(c,1) = HL_L{c,1};
end
%矸石含量可视化
plot(GS)

下面是针对R值图像，采用OTSU阈值分割，而非人工手动选择阈值，此外，采用的OTSU双阈值分割。

main.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
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
clc
clear all;
global c_im;
c_im = imread('D:\');
c_im = rgb2gray(c_im);
figure,imshow(c_im);
% c_im = c_im(1:1000,1:550);
%调用预图像处理m文件
pre;
%%
%转换数据类型并做线性缩放
R_im_new = uint8(ans);
imgB = im2double(R_im_new);
%% 先高斯滤波平滑，去高频，再用阈值处理 
gaussH=fspecial('gaussian',[5 5],10);
smoothB=imfilter(imgB,gaussH);
imshow(smoothB);
%计算双阈值
[t1B,t2B] =DoubleOtsuThresh(smoothB);    
% 如果imgB是unit8类型,则灰度范围在[0,255] 要把算法得到的两个阈值乘255
% 如果imgB是[0,1]范围的double类型则不用乘
%T1B=t1B*255;   
%T2B=t2B*255;
T1B = t1B;
T2B = t2B;
J9=img2gray(smoothB,T1B,T2B); %转为灰度图
figure; 
imshow(J9,[]);title('平滑+双阈值分割'); % imshow里面J9后面的[]相当于直方图均衡效果，让灰度值均匀分布,可以去掉对比一下
%%
%算法改进之三值图像空洞填充
% figure(12),imshow(J9,[]);
J10 = imfill(J9);
figure(13),imshow(J10,[])
J11 = J10 +J9 ;
figure(14),imshow(J11,[]);
%%
%伪色彩处理
I=double(J11);
[m,n]=size(I);
R = zeros(m,n);
G = zeros(m,n);
B = zeros(m,n);
for i=1:m
    for j=1:n
        %赋予背景层白色
               if I(i,j) ==0
                    R(i,j)=255;
                    G(i,j)=255;
                    B(i,j)=255;
               end
          %赋予中层红色
               if I(i,j)==3
                  R(i,j)=255;
                  G(i,j)=0;
                  B(i,j)=0;
               end
               if I(i,j)==2
                  R(i,j)=255;
                  G(i,j)=0;
                  B(i,j)=0;
               end
               if I(i,j)==1
                  R(i,j)=255;
                  G(i,j)=0;
                  B(i,j)=0;
               end
               %赋予矸石层绿色
               if I(i,j) == 4
                  R(i,j)=0;
                  G(i,j)=255;
                  B(i,j)=0;   
               end
           end      
end
rgbim = zeros(m,n,3);
for i=1:m
    for j=1:n
        rgbim(i,j,1)=R(i,j);
        rgbim(i,j,2)=G(i,j);
        rgbim(i,j,3)=B(i,j);
    end
end
rgbim=rgbim/256;
subplot(1,2,1);imshow(c_im);
subplot(1,2,2);
imshow(rgbim);
hold on ;
%%
%图形美化、数据处理及可视化
%标注二进制图像中已连接的部分
[L,num]=bwlabel(J9,8);     
%封装函数获取各种参数信息
stats_1 = regionprops(L,'Area','PixelList','Centroid','BoundingBox');
%绘制质心坐标
Cen_1 = cat(1,stats_1.Centroid);
figure(8),imshow(c_im);
hold on ;
for i=1:length(stats_1) 
    plot(Cen_1(i,1), Cen_1(i,2), 'r+'); 
end
%绘制最小外接矩形框
Temp_1 = cat(1,stats_1.BoundingBox);
for i=1:length(stats_1) 
    temp = Temp_1(i,:);
    rectangle('position',temp,'edgecolor','r');
end 
%标记外接矩形框起始点坐标
for i = 1:length(stats_1)
    plot(Temp_1(i,1),Temp_1(i,2),'g');
    text(stats_1(i).BoundingBox(1)-10,stats_1(i).BoundingBox(2),num2str(i),'Color','y','FontSize',12);
end
hold on;
%%
%统计各区域红色像素占比
%二值化取反图像连通域检测
[G,num] = bwlabel(J9,8);
global R_im_f;
%连通域分离
EZ = cell(num,1);
for c = 1:num
    G = bwlabel(R_im_f,8);
    for i = 1:m
        for j = 1:n
             if G(i,j) ~= c
                 G(i,j) = 0;  
             end
        end
    end
    EZ{c,1} = G .* (1/c);
end

%imshow( EZ{1,1})
HL = cell(num,1);
for c = 1:num
    HL{c,1} = J9.* EZ{c,1};
end
%figure(3),imshow(HL{18,1},[])
%计算红色像素占比
HL_L = cell(num,1);
for c = 1:num
    p = 0;
    q = 0;
    for i = 1:m
        for j = 1:n
            if  HL{c,1} (i,j) == 1
              p = p +1  ;
            end
            if HL{c,1}(i,j) == 2
                q = q+1;
            end
        end
    end
    HL_L{c,1} = p/(p+q);
end
for c= 1:num
    GS(c,1) = HL_L{c,1};
end


%%视情况而变
%可视化
figure(6),plot(GS(1:num,1),'r+');
hold on;

plot(GS(31:61,1),'g+'),hold on;

plot(GS(1:30,1));
hold on;
plot(GS(31:61,1))

figure(7),scatter(1:num,GS,'r+');

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
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
clc
clear all;
%%  图B 彩色图
load('R_im_new.mat');
%返回数据类型double 范围【0,1】
imgB = im2double(R_im_new);
%R_im_new = bwareaopen(R_im_new,10);
%% 先高斯滤波平滑图B，去高频，再用阈值处理 
gaussH=fspecial('gaussian',[5 5],10);
smoothB=imfilter(imgB,gaussH);
[t1B,t2B] =DoubleOtsuThresh(smoothB);    
% 如果imgB是unit8类型,则灰度范围在[0,255] 要把算法得到的两个阈值乘255
% 如果imgB是[0,1]范围的double类型则不用乘
%T1B=t1B*255;   
%T2B=t2B*255;
T1B=t1B;  
T2B=t2B;
J9=img2gray(smoothB,0.1725,0.5020); %转为灰度图
figure 
imshow(J9,[]);title('平滑+双阈值分割'); % imshow里面J9后面的[]相当于直方图均衡效果，让灰度值均匀分布,可以去掉对比一下
%%
%伪色彩处理
I=double(J9);
[m,n]=size(I);
L=256;
for i=1:m
    for j=1:n
        if I(i,j)==0
             R(i,j)=255;
             G(i,j)=255;
              B(i,j)=255;
        else if I(i,j)==1
             R(i,j)=0;
                G(i,j)=255;
                B(i,j)=0;
        else if I(i,j)==2
            R(i,j)=255;
            G(i,j)=0;
            B(i,j)=0;
            end
            end
        end
    end
end

for i=1:m
    for j=1:n
        rgbim(i,j,1)=R(i,j);
        rgbim(i,j,2)=G(i,j);
        rgbim(i,j,3)=B(i,j);
    end
end

rgbim=rgbim/256;

subplot(1,2,2);
imshow(rgbim);

%%
%图形美化
%标注二进制图像中已连接的部分
[L,num]=bwlabel(Low_im_bw,8);     
%封装函数获取各种参数信息
stats_1 = regionprops(L,'Area','PixelList','Centroid','BoundingBox');
%绘制质心坐标
Cen_1 = cat(1,stats_1.Centroid);
subplot(1,2,1);imshow(Low_im);
hold on ;
for i=1:length(stats_1) 
    plot(Cen_1(i,1), Cen_1(i,2), 'r+'); 
end
%绘制最小外接矩形框
Temp_1 = cat(1,stats_1.BoundingBox);
for i=1:length(stats_1) 
    temp = Temp_1(i,:);
    rectangle('position',temp,'edgecolor','r');
end 
%标记外接矩形框起始点坐标
for i = 1:length(stats_1)
    plot(Temp_1(i,1),Temp_1(i,2),'c');
    text(stats_1(i).BoundingBox(1)-10,stats_1(i).BoundingBox(2),num2str(i),'Color','y','FontSize',12);
end
hold on;

img2gray.m:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
function out=img2gray(img,t1,t2)
 [M,N]=size(img);
 out=zeros(M,N);
  
 for m=1:M
     for n=1:N
         if img(m,n)<t1
             out(m,n)=0;
         else
             if( img(m,n)>t1 && img(m,n)<t2)
                 out(m,n)=1;  %灰度级可调
             else
                out(m,n)=2;	%灰度级可调
             end
         end
     end
 end
end

调用函数：DoubleOtsuThresh.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
45
46
47
48
49
50
51
function [t1,t2]=DoubleOtsuThresh(img)
%#####%
%  Otsu 双阈值求解
%  输入 图像img，输出最优阈值t1和t2(归一化,范围在[0,1])
%#####%
BinsNum = 256;  
hist = imhist(img,BinsNum);   
p = hist / sum(hist);          % 直方图的概率密度函数
mG= sum(p .* (1:BinsNum)');     % 全局均值

P1 = cumsum(p);                 % 概率分布
m1 = cumsum(p .* (1:BinsNum)')./P1; % 256*1  每个阈值的前景平均灰度

% 根据算法理论，从k2+1累加到L-1，可以先倒着累加再翻转回来
P3= cumsum(flip(p));          
m3 = cumsum(flip(p) .* flip(1:BinsNum)')./P3;
P3=flip(P3);
P3=[P3(2:end) ;0];             %P3的索引用k2，则P3(1)实际上是从p(2)加到p(256)，所以要移动一个
%移动后也用不到P3(1)这个值，因为k2>k1>1
m3=flip(m3);
m3=[m3(2:end) ;0];

% m2 P2为 k1*k2的索引矩阵，为 256*256 的上三角矩阵 因为k1<k2
m2=zeros(BinsNum,BinsNum);    
P2=zeros(BinsNum,BinsNum);

for k1=1:BinsNum-2
    for k2=k1+1:BinsNum-1
        P2(k1,k2)=1-P1(k1)-P3(k2);
        m2(k1,k2) = sum( (k1+1:k2)' .* p(k1+1:k2) )./ P2(k1,k2) ;        
    end
end

% 遍历k1,k2各种组合求方差
variance=zeros(BinsNum,BinsNum);
for k1=1:BinsNum-2
    for k2=k1+1:BinsNum-1
        
% variance 为256*256 的上三角矩阵 因为l1<k2
    variance(k1,k2) = P1(k1)*(m1(k1)-mG)^2 + ...
                      P2(k1,k2)*(m2(k1,k2)-mG)^2 + ...
                      P3(k2)*(m3(k2)-mG)^2;
    end
end
% 最大方差点即为最优阈值
[~,index] = max(variance(:));
[index_row,index_col] = ind2sub(size(variance),index);
%灰度值为[0,255] 因此需要-1
t1= (index_row-1)./(BinsNum-1);
t2= (index_col-1)./(BinsNum-1);
end

采用双阈值分割，并赋予图像三种颜色：doubleTH.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
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
clc 
clear all;
%读取高低能图像
% High_im = imread('D:\双能X射线\煤(20-50mm)\2020_7_31_9_22_13_928_H.bmp');
% Low_im = imread('D:\双能X射线\煤(20-50mm)\2020_7_31_9_22_13_928_L.bmp');
High_im = imread('D:\双能X射线\50-180mm厚度煤矸混合\2020_8_5_14_34_22_760_H.bmp');
Low_im = imread('D:\双能X射线\50-180mm厚度煤矸混合\2020_8_5_14_34_22_760_L.bmp');
% High_im = imread('D:\双能X射线\矸石(20-50mm)\2020_7_31_9_36_46_732_H.bmp');
% Low_im = imread('D:\双能X射线\矸石(20-50mm)\2020_7_31_9_36_46_732_L.bmp');

High_im = High_im(40:1000,80:1000);
Low_im = Low_im(40:1000,80:1000);
figure,
subplot(1,2,1),imshow(High_im);
subplot(1,2,2),imshow(Low_im);
%中值滤波
High_im = medfilt2(High_im,[3,3]);
Low_im = medfilt2(Low_im,[3,3]);
%计算二值化图像
thersh_H = graythresh(High_im);
thersh_L = graythresh(Low_im);
High_im_bw = im2bw(High_im,thersh_H);
Low_im_bw = im2bw(Low_im,thersh_L);
%二值化取反
High_im_bw = ~High_im_bw;
Low_im_bw = ~Low_im_bw ;
% High_im_bw = bwareaopen(High_im_bw,5);
% Low_im_bw = bwareaopen(Low_im_bw,5);
figure,
subplot(1,2,1),imshow(High_im_bw);
subplot(1,2,2),imshow(Low_im_bw);

% Low_im_nobord = imclearborder(Low_im_bw,4);
% seD = strel('diamond',1);
% BWfinal = imerode(Low_im_nobord,seD);
% figure,
%  imshow(BWfinal);
%空皮带平均灰度值构建矩阵
% em_im_H = imread('D:\2020_7_31_9_30_26_154_H.bmp');
% em_im_L= imread('D:\2020_7_31_9_30_26_154_L.bmp');
% em_im_H = em_im_H(300:300,400:400);
% mean_H = mean(mean(em_im_H));
% em_im_L = em_im_L(300:300,400:400);
% mean_L = mean(mean(em_im_L));
em_H = double(High_im.*(0)) +191;
em_L = double(Low_im.*(0)) +195;
%计算并显示R值图像
R_im= (log(double(em_L) ./ double(Low_im)) )./ (log(double(em_H) ./ double(High_im)));
figure,
imshow(R_im);
%掩膜处理提取目标区域R值图像
R_im_f = R_im .* Low_im_bw;
figure,
imshow(R_im_f);

R_im_f(find(isnan(R_im_f)==1)) = 0;
%%
%映射R值图像到0~255灰度区间内
Omax = 255;
Omin = 0;
x = find(R_im_f <1 & R_im_f >0.8 );
y = find(R_im_f <1.6 & R_im_f >0);
p = R_im_f (x);
q = R_im_f (y);
a =min(p);
b =max(q);
I = ones(961,921);
for i = 1 : 961
    for j = 1 : 921
        R_im_new(i,j) =(Omax - Omin) * (R_im_f(i,j) - a) / (b - a) + Omin;
        if R_im_new(i,j) == 0
            R_im_new(i,j) =255;
        end
%         if R_im_new(i,j) > 255
%             R_im_new(i,j) =0;
%         end
    end
end      
figure(6),imshow(uint8(R_im_new),[0,255]);

figure,imagesc(uint8(R_im_new),[0,255]);colorbar
colormap jet ;
%%

% R_im_f = uint8(R_im_f);
% figure,imhist(R_im_f);
%将R值图像中NaN置零
R_im_f(find(isnan(R_im_f)==1)) = 0;

%阈值分割
img = R_im_f;
[M,N]=size(img);
out=zeros(M,N);
 
t1 =1.25;   
t2 = 1.3;  

 for m=1:M
     for n=1:N
         if (img(m,n)>=t1 && img(m,n)<=t2 )
             out(m,n)=3;
         end
         if(img(m,n)<t1)
                 out(m,n)=1;  %灰度级可调
         end
         if(img(m,n) == 0)
                 out(m,n)=0;  %灰度级可调
         end
         if(img(m,n) > t2)
                 out(m,n)=2;  %灰度级可调
         end
      end
 end
figure,imshow(out);
 %伪色彩处理
I=double(out);
[m,n]=size(I);
R = zeros(m,n);
G = zeros(m,n);
B = zeros(m,n);
for i=1:m
    for j=1:n
        %赋予背景层白色
               if I(i,j) ==0
                    R(i,j)=255;
                    G(i,j)=255;
                    B(i,j)=255;
               end
         %赋予煤层绿色
               if I(i,j)==1
                  R(i,j)=0;
                  G(i,j)=255;
                  B(i,j)=0;
               end
         %赋予矸石层红色             
               if I(i,j)==2
                  R(i,j)=255;
                  G(i,j)=0;
                  B(i,j)=0;
               end
          %赋予中间层蓝色             
               if I(i,j)==3
                  R(i,j)=0;
                  G(i,j)=0;
                  B(i,j)=255;
               end
           end      
end

rgbim = zeros(m,n,3);
for i=1:m
    for j=1:n
        rgbim(i,j,1)=R(i,j);
        rgbim(i,j,2)=G(i,j);
        rgbim(i,j,3)=B(i,j);
    end
end

rgbim=rgbim / 256;
subplot(1,2,1);imshow(High_im);
subplot(1,2,2);
imshow(rgbim);
hold on ;

%统计各区域红色像素占比
%二值化取反图像连通域检测
[G,num] = bwlabel(Low_im_bw);
%连通域分离
EZ = cell(num,1);
for c = 1:num
    G = bwlabel(Low_im_bw,8);
    for i = 1:m
        for j = 1:n
             if G(i,j) ~= c
                 G(i,j) = 0;  
             end
        end
    end
    EZ{c,1} = G .* (1/c);
end
%imshow( EZ{1,1})
HL = cell(num,1);

for c = 1:num
    HL{c,1} = out.* EZ{c,1};
end
%figure,imshow(HL{21,1},[])
%计算红色像素占比
HL_L = cell(num,1);
for c = 1:num
    p = 0;
    q = 0;
    for i = 1:m
        for j = 1:n
            if  HL{c,1} (i,j) == 1
              p = p +1  ;
            end
            if HL{c,1}(i,j) == 2
                q = q+1;
            end
        end
    end
    HL_L{c,1} = q/(p+q);
end
for c= 1:num
    GS(c,1) = HL_L{c,1};
end
%矸石含量可视化
plot(GS)

R值阈值图像分割结果，处理后的图像：

Snipaste_2024-05-08_10-24-41.png