利用Matlab仿真实现图像烟雾识别(k-means聚类图像分割+LBP+PCA+SVM)

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
%基于LBP特征提取的主程序代码
clc; 
clear ;  
k = 2;
acc1 = 0;
acc2 = 0;
acc = 0;
%%  标签制作  
ReadList1  = textread('pos_list.txt','%s','delimiter','n');%载入正样本列表  
sz1=size(ReadList1);   
label1=ones(sz1(1),1); %正样本标签  
ReadList2  = textread('neg_list.txt','%s','delimiter','n');%载入负样本列表
sz2=size(ReadList2);  
label2=zeros(sz2(1),1);%负样本标签  
label_train = [label1',label2'];%训练集标签
ReadList_pos = textread('pos_test_list.txt','%s','delimiter','n');%载入测试正样本列表  
sz_pos=size(ReadList_pos);   
label_pos=ones(sz_pos(1),1); %正样本标签
ReadList_neg  = textread('neg_test_list.txt','%s','delimiter','n');%载入测试负样本列表
sz_neg=size(ReadList_neg);  
label_neg=zeros(sz_neg(1),1);%负样本标签  
label_test = [label_pos',label_neg'];%测试集误差
total_trainnum=length(label_train);  
total_testnum = length(label_test);
data1 = zeros(total_trainnum,256);  
data2 = zeros(total_testnum,256);
%% 提取特征
%读取训练集正样本并计算lbp特征 
for i=1:sz1(1)
   name=char(ReadList1(i,1));  
   image1=imread(strcat('F:模式识别matlab程序模式识别大作业yanwujiancepos',name));
    I=double(image1)/255;
   clu_kmeans=imkmeans(I,3);
   clu_pic=clu_kmeans/3;
   lbps = lbp(clu_pic);
   data1(i,:)=lbps;  
end
%读取训练集负样本并计算lbp特征  
for j=1:sz2(1)
   name= char(ReadList2(j,1));  
   image2=imread(strcat('F:模式识别matlab程序模式识别大作业yanwujianceneg',name));  
    I=double(image2)/255;
   clu_kmeans=imkmeans(I,3);
   clu_pic=clu_kmeans/3;
   lbps = lbp(clu_pic);
   data1(sz1(1)+j,:)=lbps;  
end
%读取测试集正样本并计算lbp特征
for m=1:sz_pos(1)
   test_name= char(ReadList_pos(m,1));  
   image3=imread(strcat('F:模式识别matlab程序模式识别大作业yanwujiancetestpos_test',test_name));  
    I=double(image3)/255;
   clu_kmeans=imkmeans(I,3);
   clu_pic=clu_kmeans/3;
   lbpst= lbp(clu_pic);
   data2(m,:)=lbpst;  
end
%读取测试集负样本并计算lbp特征
for n =1:sz_neg(1)
    test_name=char(ReadList_neg(n,1)); 
    image4=imread(strcat('F:模式识别matlab程序模式识别大作业yanwujiancetestneg_test',test_name));
     I=double(image4)/255;
   clu_kmeans=imkmeans(I,3);
   clu_pic=clu_kmeans/3;
   lbps = lbp(clu_pic);
    data2(sz_pos(1)+n,:)=lbpst; 
end
load data1
load data2
load svmStruct3
%数据降维
[COEFF SCORE latent]=princomp(data1(:,:));%训练集数据降维
pcaData1 = SCORE(:,1:k);
latent = 100*latent/sum(latent);
for i = 1:8
latent(i+1) = latent(i+1)+latent(i)
end
plot(latent(1:8));%画出前8个特征值所包含的图像信息比例
x0 = bsxfun(@minus,data2,mean(data2,1));
pcaData2_sw = x0*COEFF(:,:);
pcaData2 = pcaData2_sw(:,1:k);
%%  评估方法：交叉验证法
[train, test] = crossvalind('holdOut',label_train);   %随机选择训练集合测试集
cp = classperf(label_train);  %评估分类器性能
svmStruct3hog = svmtrain(pcaData1(train,1:k),label_train(train));%训练SVM分类器  
%使用svmtrain进行训练,得到训练后的结构svmStruct3hog,在预测时使用
save svmStruct3hog   %%保存 svmStruct3hog
cros = svmclassify(svmStruct3hog,pcaData1(test,1:k)); 
classperf(cp,cros ,test);  
cp.CorrectRate   
%% 测试
load svmStruct3hog
for i=1:sz_pos(1)
       classes = svmclassify(svmStruct3,pcaData2(i,:));%classes的值即为分类结果
       if classes==1
           acc1=acc1+1;%记录正确分类的样本数
       end
end
for j = sz_pos(1)+1:1383
       classes = svmclassify(svmStruct3,pcaData2(j,:));%classes的值即为分类结果
       if classes~=1
           acc2=acc2+1;%记录正确分类的样本数
       end
end 
acc = acc1+acc2;
fprintf('精确度为：%5.2f%%n',(acc/(sz_neg(1)+sz_pos(1)))*100);%计算预测的正确率

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
%lbp特征提取代码
function result = lbp(varargin) % image,radius,neighbors,mapping,mode)
% Check number of input arguments.
error(nargchk(1,5,nargin));
image=varargin{1};
d_image=double(image);
 
if nargin==1
    spoints=[-1 -1; -1 0; -1 1; 0 -1; -0 1; 1 -1; 1 0; 1 1];
    neighbors=8;
    mapping=0;
    mode='h';
end
 
if (nargin == 2) && (length(varargin{2}) == 1)
    error('Input arguments');
end
 
if (nargin > 2) && (length(varargin{2}) == 1)
    radius=varargin{2};
    neighbors=varargin{3};
    spoints=zeros(neighbors,2);
 
    % Angle step.
    a = 2*pi/neighbors;
    for i = 1:neighbors
        spoints(i,1) = -radius*sin((i-1)*a);
        spoints(i,2) = radius*cos((i-1)*a);
    end
   
    if(nargin >= 4)
        mapping=varargin{4};
        if(isstruct(mapping) && mapping.samples ~= neighbors)
            error('Incompatible mapping');
        end
    else
        mapping=0;
    end
   
    if(nargin >= 5)
        mode=varargin{5};
    else
        mode='h';
    end
end
 
 
if (nargin > 1) && (length(varargin{2}) > 1)
    spoints=varargin{2};
    neighbors=size(spoints,1);
   
    if(nargin >= 3)
        mapping=varargin{3};
        if(isstruct(mapping) && mapping.samples ~= neighbors)
            error('Incompatible mapping');
        end
    else
        mapping=0;
    end
   
    if(nargin >= 4)
        mode=varargin{4};
    else
        mode='h';
    end  
end
 
% Determine the dimensions of the input image.
[ysize xsize] = size(image);
 
miny=min(spoints(:,1));
maxy=max(spoints(:,1));
minx=min(spoints(:,2));
maxx=max(spoints(:,2));
 
% Block size, each LBP code is computed within a block of size bsizey*bsizex
bsizey=ceil(max(maxy,0))-floor(min(miny,0))+1;
bsizex=ceil(max(maxx,0))-floor(min(minx,0))+1;
 
 
% Coordinates of origin (0,0) in the block
origy=1-floor(min(miny,0));
origx=1-floor(min(minx,0));
 
 
% Minimum allowed size for the input image depends
% on the radius of the used LBP operator.
if(xsize < bsizex || ysize < bsizey)
  error('Too small input image. Should be at least (2*radius+1) x (2*radius+1)');
end
 
 
% Calculate dx and dy;
dx = xsize - bsizex;
dy = ysize - bsizey;
 
 
% Fill the center pixel matrix C.
C = image(origy:origy+dy,origx:origx+dx);
d_C = double(C);
 
 
bins = 2^neighbors;
 
 
% Initialize the result matrix with zeros.
result=zeros(dy+1,dx+1);
 
 
%Compute the LBP code image
 
 
for i = 1:neighbors
  y = spoints(i,1)+origy;
  x = spoints(i,2)+origx;
  % Calculate floors, ceils and rounds for the x and y.
  fy = floor(y); cy = ceil(y); ry = round(y);
  fx = floor(x); cx = ceil(x); rx = round(x);
  % Check if interpolation is needed.
  if (abs(x - rx) < 1e-6) && (abs(y - ry) < 1e-6)
    % Interpolation is not needed, use original datatypes
    N = image(ry:ry+dy,rx:rx+dx);
    D = N >= C;
  else
    % Interpolation needed, use double type images
    ty = y - fy;
    tx = x - fx;
 
 
    % Calculate the interpolation weights.
    w1 = (1 - tx) * (1 - ty);
    w2 =      tx  * (1 - ty);
    w3 = (1 - tx) *      ty ;
    w4 =      tx  *      ty ;
    % Compute interpolated pixel values
    N = w1*d_image(fy:fy+dy,fx:fx+dx) + w2*d_image(fy:fy+dy,cx:cx+dx) + ...
        w3*d_image(cy:cy+dy,fx:fx+dx) + w4*d_image(cy:cy+dy,cx:cx+dx);
    D = N >= d_C;
  end 
  % Update the result matrix.
  v = 2^(i-1);
  result = result + v*D;
end
 
 
%Apply mapping if it is defined
if isstruct(mapping)
    bins = mapping.num;
    for i = 1:size(result,1)
        for j = 1:size(result,2)
            result(i,j) = mapping.table(result(i,j)+1);
        end
    end
end
 
 
if (strcmp(mode,'h') || strcmp(mode,'hist') || strcmp(mode,'nh'))
    % Return with LBP histogram if mode equals 'hist'.
    result=hist(result(:),0:(bins-1));
    if (strcmp(mode,'nh'))
        result=result/sum(result);
    end
else
    %Otherwise return a matrix of unsigned integers
    if ((bins-1)<=intmax('uint8'))
        result=uint8(result);
    elseif ((bins-1)<=intmax('uint16'))
        result=uint16(result);
    else
        result=uint32(result);
    end
end
end

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
%k-means图像聚类分割
function [F,C]=imkmeans(I,C)
% I:图像矩阵,支持彩色或者灰度图
% C:聚类中心，可以是整数或者数组，整数表示随机选择K个聚类中心
% F:样本聚类编号
if nargin~=2
    error('IMKMEANS:InputParamterNotRight','只能有两个输入参数！');
end
if isempty(C)
    K=2;
    C=[];
elseif isscalar(C)
    K=C;
    C=[];
else
    K=size(C,1);
end
%% I.提取像素点特征向量
X=exactvecotr(I);
%% II.搜索初始聚类中心
if isempty(C)
    C=searchintial(X,'sample',K);
end
%% III.循环搜索聚类中心
Cprev=rand(size(C));
while true
    %计算样本到中心的距离
    D=sampledist(X,C,'euclidean');
    %找出最近的聚类中心
    [~,locs]=min(D,[],2);
    %使用样本均值更新中心
    for i=1:K
        C(i,:)=mean(X(locs==i,:),1);
    end
    %判断聚类算法是否收敛
    if norm(C(:)-Cprev(:))<eps
        break
    end
    %保存上一次聚类中心
    Cprev=C;
end
[m,n,~]=size(I);
F=reshape(locs,[m,n]);