一、算法簡(jiǎn)介
1.1 c-means聚類算法
聚類分析是根據(jù)在數(shù)據(jù)中發(fā)現(xiàn)的描述對(duì)象及其關(guān)系的信息，將數(shù)據(jù)對(duì)象進(jìn)行分組。目的是使組內(nèi)的對(duì)象相互之間是相似的（相關(guān)的），而不同組中的對(duì)象是不同的（不相關(guān)的）。組內(nèi)相似性越大，組間差距越大，說明聚類效果越好。
也就是說，聚類的目標(biāo)是得到較高的類內(nèi)相似度和較低的類間相似度，使得類間的距離盡可能大，類內(nèi)樣本與類中心的距離盡可能小。在此，我們選用k-means聚類算法。
1 .2 LBP算法
LBP（Local Binary Pattern，局部二值模式）是一種用來描述圖像局部紋理特征的算子；它具有旋轉(zhuǎn)不變性和灰度不變性等顯著的優(yōu)點(diǎn)。它是首先由T. Ojala, M.Pietikäinen, 和D. Harwood 在1994年提出，用于紋理特征提取，提取的特征是圖像的局部的紋理特征。
原始的LBP算子定義為在3*3的窗口內(nèi)，以窗口中心像素為閾值，將相鄰的8個(gè)像素的灰度值與其進(jìn)行比較，若周圍像素值大于中心像素值，則該像素點(diǎn)的位置被標(biāo)記為1，否則為0。這樣，3*3鄰域內(nèi)的8個(gè)點(diǎn)經(jīng)比較可產(chǎn)生8位二進(jìn)制數(shù)（通常轉(zhuǎn)換為十進(jìn)制數(shù)即LBP碼，共256種），即得到該窗口中心像素點(diǎn)的LBP值，并用這個(gè)值來反映該區(qū)域的紋理信息。
1.3 PCA算法
PCA(Principal Component Analysis)，即主成分分析方法，是一種使用最廣泛的數(shù)據(jù)降維算法。其算法步驟如下：
1）數(shù)據(jù)中心化——去均值，根據(jù)需要，有的需要?dú)w一化——Normalized；
2）求解協(xié)方差矩陣；
3）利用特征值分解/奇異值分解 求解特征值以及特征向量；
4）將特征值從大到小排序，保留前k個(gè)特征向量
5）利用特征向量構(gòu)造投影矩陣；
6）利用投影矩陣，得出降維的數(shù)據(jù)。
1.4 SVM算法
支持向量機(jī)（support vector machines, SVM）是一種二分類模型，它的基本模型是定義在特征空間上的間隔最大的線性分類器，間隔最大使它有別于感知機(jī)；SVM還包括核技巧，這使它成為實(shí)質(zhì)上的非線性分類器。SVM的的學(xué)習(xí)策略就是間隔最大化，可形式化為一個(gè)求解凸二次規(guī)劃的問題，也等價(jià)于正則化的合頁損失函數(shù)的最小化問題。SVM的的學(xué)習(xí)算法就是求解凸二次規(guī)劃的最優(yōu)化算法。
SVM學(xué)習(xí)的基本想法是求解能夠正確劃分訓(xùn)練數(shù)據(jù)集并且?guī)缀伍g隔最大的分離超平面。如下圖所示即為分類超平面，對(duì)于線性可分的數(shù)據(jù)集來說，這樣的超平面有無窮多個(gè)（即感知機(jī)），但是幾何間隔最大的分類超平面卻是唯一的。如下圖1-1SVM算法示意圖
圖1-1SVM算法示意圖
二、算法實(shí)現(xiàn)
2.1 煙霧識(shí)別算法流程
1）首先對(duì)所有圖像進(jìn)行預(yù)處理，假定將有煙當(dāng)作正樣本，將沒煙看作負(fù)樣本，train集的smoke文件夾改名為pos，train集的non文件夾改名為neg；同理將test集的smoke文件夾改名為pos，test集的non文件夾改名為neg。為了對(duì)所有圖片進(jìn)行處理，將train和test中的pos和neg中的圖片全部規(guī)范命名格式為0001.jpg、0002.jpg、0003.jpg、0004.jpg、0005.jpg......。將這些圖片名字提取出來分別存到“pos_list.txt、neg_list.txt、pos_test_list.txt、neg_test_list.txt文本中。如下圖2-1圖2-2所示
圖2-1
圖2-2
2）利用c-means聚類算法對(duì)訓(xùn)練集和測(cè)試集圖像的像素進(jìn)行聚類，實(shí)現(xiàn)圖像分割。
3）利用LBP對(duì)分割后的訓(xùn)練集圖像和測(cè)試集圖像進(jìn)行特征提取。
4）分別對(duì)訓(xùn)練集和測(cè)試集使用主成分分析法（PCA）進(jìn)行特征降維。
5）利用對(duì)訓(xùn)練集降維后得到的二維特征訓(xùn)練SVM二分類模型，
6）最后利用對(duì)測(cè)試集降維后得到的二維特征進(jìn)行分類預(yù)測(cè)。
整體算法流程如下圖2-3所示
圖2-3 算法流程框圖
2.2 c-means算法實(shí)現(xiàn)
圖像分割是利用圖像的灰度、顏色、紋理、形狀等特征，把圖像分成若干個(gè)互不重疊的區(qū)域，并使這些特征在同一區(qū)域內(nèi)呈現(xiàn)相似性，在不同的區(qū)域之間存在明顯的差異性。然后就可以將分割的圖像中具有獨(dú)特性質(zhì)的區(qū)域提取出來用于不同的研究。圖像識(shí)別的基礎(chǔ)是圖像分割，其作用是把反映物體真實(shí)情況的、占據(jù)不同區(qū)域的、具有不同特性的目標(biāo)區(qū)分開來，并形成數(shù)字特征。因此本文利用c-means聚類算法實(shí)現(xiàn)圖像分割，實(shí)現(xiàn)對(duì)噪聲的過濾，在構(gòu)建煙霧識(shí)別模型的過程中，首先分別對(duì)無煙和有煙的圖像進(jìn)行c-means聚類圖像分割。
本文對(duì)預(yù)處理過后的訓(xùn)練集和測(cè)試集圖像進(jìn)行像素聚類，在此分別列舉一張有煙圖和無煙圖的圖像分割前后的效果對(duì)比。如圖2-4和圖2-5所示
圖2-4 無煙圖像分割前后對(duì)照?qǐng)D
圖2-5有煙圖像分割前后對(duì)照?qǐng)D
2.3 LBP算法實(shí)現(xiàn)
本文LBP算法將像素聚類（3類）以后的圖像進(jìn)行特征提取。在此分別列舉一張有煙圖和無煙圖的圖像特征提取前后的效果對(duì)比。
圖2-6無煙三像素聚類LBP特征提取前后對(duì)照?qǐng)D
圖2-7有煙三像素聚類LBP特征提取前后對(duì)照?qǐng)D
本文PCA算法將HOG或LBP提取的特征進(jìn)行特征降維，使數(shù)據(jù)可視化。PCA算法可以獲取原有特征的大部分信息，降維以后的前k個(gè)特征值保留下來的信息占原有信息的比例可有下式計(jì)算獲得。
對(duì)LBP算法提取的特征進(jìn)行特征降維，在此取前兩維特征進(jìn)行模型訓(xùn)練，前兩維度保留的信息含有98.75%，如下圖2-8所示.
2.4 SVM算法實(shí)現(xiàn)
在經(jīng)過上述圖像預(yù)處理、圖像像素聚類、LBP特征提取、PCA特征降維至兩維過程之后，將二維特征向量作為輸入訓(xùn)練SVM模型，最終得到模型在訓(xùn)練集上的分類準(zhǔn)確度。
利用k-means+LBP+PCA+SVM算法，多次訓(xùn)練模型，最終取平均值，得到在訓(xùn)練集上的分類準(zhǔn)確度為79%，在測(cè)試集上的分類準(zhǔn)確度為78%。下圖為模型在訓(xùn)練集上的分類效果圖。
三、結(jié)果分析
經(jīng)過第二章的算法實(shí)現(xiàn)，最終得到了完整的SVM二分類模型，利用該模型對(duì)test中的pos樣本的圖片和neg樣本的圖片進(jìn)行預(yù)測(cè)。預(yù)測(cè)前，首先需要對(duì)測(cè)試集圖片經(jīng)過預(yù)處理、其次利用k-means3聚類法對(duì)像素進(jìn)行聚類得到最終圖像分割聚類圖、然后對(duì)聚類圖進(jìn)行LBP特征提取、最后再利用PCA對(duì)提取出來的特征進(jìn)行特征降維。將最終得到的二維特征向量作為模型的輸入，進(jìn)行分類預(yù)測(cè)，最終得到結(jié)果。對(duì)于LBP特征提取方法，在訓(xùn)練集和測(cè)試集上的準(zhǔn)確率分別為79%和78%。經(jīng)過對(duì)比可以發(fā)現(xiàn)模型的泛化性能良好。
最后筆者不得不提的是，之所以采取上訴方法實(shí)現(xiàn)煙霧識(shí)別是因?yàn)椋笞鳂I(yè)要求必須包含聚類、分類、降維。筆者也嘗試過直接使用LBP+SVM實(shí)現(xiàn)煙霧識(shí)別的方法，并且對(duì)測(cè)試集的準(zhǔn)確率可以達(dá)到93%。
這是兩種不一樣的解決問題的思路。若采用本文的思路是Pipeline，若直接采用LBP+SVM的思路叫做end2end，各有優(yōu)缺點(diǎn)。Pipeline是將一個(gè)問題拆解成若干個(gè)子問題一次解決，然后串在一起，這種方法易于實(shí)現(xiàn)，且靈活性和可解釋性更高，但缺點(diǎn)是多個(gè)子任務(wù)會(huì)造成錯(cuò)誤累積。end2end是將一個(gè)問題看成一個(gè)整體，一般可以獲得比pipeline更高的性能，但是整體像一個(gè)黑盒，可解釋性差。現(xiàn)在深度學(xué)習(xí)最新研究的趨勢(shì)是end2end的方法。
%基于LBP特征提取的主程序代碼
clc;
clear ;
k = 2;
acc1 = 0;
acc2 = 0;
acc = 0;
%% 標(biāo)簽制作
ReadList1 = textread('pos_list.txt','%s','delimiter','\n');%載入正樣本列表
sz1=size(ReadList1);
label1=ones(sz1(1),1); %正樣本標(biāo)簽
ReadList2 = textread('neg_list.txt','%s','delimiter','\n');%載入負(fù)樣本列表
sz2=size(ReadList2);
label2=zeros(sz2(1),1);%負(fù)樣本標(biāo)簽
label_train = [label1',label2'];%訓(xùn)練集標(biāo)簽
ReadList_pos = textread('pos_test_list.txt','%s','delimiter','\n');%載入測(cè)試正樣本列表
sz_pos=size(ReadList_pos);
label_pos=ones(sz_pos(1),1); %正樣本標(biāo)簽
ReadList_neg = textread('neg_test_list.txt','%s','delimiter','\n');%載入測(cè)試負(fù)樣本列表
sz_neg=size(ReadList_neg);
label_neg=zeros(sz_neg(1),1);%負(fù)樣本標(biāo)簽
label_test = [label_pos',label_neg'];%測(cè)試集誤差
total_trainnum=length(label_train);
total_testnum = length(label_test);
data1 = zeros(total_trainnum,256);
data2 = zeros(total_testnum,256);
%% 提取特征
%讀取訓(xùn)練集正樣本并計(jì)算lbp特征
for i=1:sz1(1)
name=char(ReadList1(i,1));
image1=imread(strcat('F:\模式識(shí)別matlab程序\模式識(shí)別大作業(yè)\yanwujiance\pos\',name));
I=double(image1)/255;
clu_kmeans=imkmeans(I,3);
clu_pic=clu_kmeans/3;
lbps = lbp(clu_pic);
data1(i,:)=lbps;
end
%讀取訓(xùn)練集負(fù)樣本并計(jì)算lbp特征
for j=1:sz2(1)
name= char(ReadList2(j,1));
image2=imread(strcat('F:\模式識(shí)別matlab程序\模式識(shí)別大作業(yè)\yanwujiance\neg\',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
%讀取測(cè)試集正樣本并計(jì)算lbp特征
for m=1:sz_pos(1)
test_name= char(ReadList_pos(m,1));
image3=imread(strcat('F:\模式識(shí)別matlab程序\模式識(shí)別大作業(yè)\yanwujiance\test\pos_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
%讀取測(cè)試集負(fù)樣本并計(jì)算lbp特征
for n =1:sz_neg(1)
test_name=char(ReadList_neg(n,1));
image4=imread(strcat('F:\模式識(shí)別matlab程序\模式識(shí)別大作業(yè)\yanwujiance\test\neg_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
%數(shù)據(jù)降維
[COEFF SCORE latent]=princomp(data1(:,:));%訓(xùn)練集數(shù)據(jù)降維
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個(gè)特征值所包含的圖像信息比例
x0 = bsxfun(@minus,data2,mean(data2,1));
pcaData2_sw = x0*COEFF(:,:);
pcaData2 = pcaData2_sw(:,1:k);
%% 評(píng)估方法：交叉驗(yàn)證法
[train, test] = crossvalind('holdOut',label_train); %隨機(jī)選擇訓(xùn)練集合測(cè)試集
cp = classperf(label_train); %評(píng)估分類器性能
svmStruct3hog = svmtrain(pcaData1(train,1:k),label_train(train));%訓(xùn)練SVM分類器
%使用svmtrain進(jìn)行訓(xùn)練,得到訓(xùn)練后的結(jié)構(gòu)svmStruct3hog,在預(yù)測(cè)時(shí)使用
save svmStruct3hog %%保存 svmStruct3hog
cros = svmclassify(svmStruct3hog,pcaData1(test,1:k));
classperf(cp,cros ,test);
cp.CorrectRate
%% 測(cè)試
load svmStruct3hog
for i=1:sz_pos(1)
classes = svmclassify(svmStruct3,pcaData2(i,:));%classes的值即為分類結(jié)果
if classes==1
acc1=acc1+1;%記錄正確分類的樣本數(shù)
end
end
for j = sz_pos(1)+1:1383
classes = svmclassify(svmStruct3,pcaData2(j,:));%classes的值即為分類結(jié)果
if classes~=1
acc2=acc2+1;%記錄正確分類的樣本數(shù)
end
end
acc = acc1+acc2;
fprintf('精確度為：%5.2f%%\n',(acc/(sz_neg(1)+sz_pos(1)))*100);%計(jì)算預(yù)測(cè)的正確率
%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
%k-means圖像聚類分割
function [F,C]=imkmeans(I,C)
% I:圖像矩陣,支持彩色或者灰度圖
% C:聚類中心，可以是整數(shù)或者數(shù)組，整數(shù)表示隨機(jī)選擇K個(gè)聚類中心
% F:樣本聚類編號(hào)
if nargin~=2
error('IMKMEANS:InputParamterNotRight','只能有兩個(gè)輸入?yún)?shù)！');
end
if isempty(C)
K=2;
C=[];
elseif isscalar(C)
K=C;
C=[];
else
K=size(C,1);
end
%% I.提取像素點(diǎn)特征向量
X=exactvecotr(I);
%% II.搜索初始聚類中心
if isempty(C)
C=searchintial(X,'sample',K);
end
%% III.循環(huán)搜索聚類中心
Cprev=rand(size(C));
while true
%計(jì)算樣本到中心的距離
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]);
以上就是利用Matlab仿真實(shí)現(xiàn)圖像煙霧識(shí)別(k-means聚類圖像分割+LBP+PCA+SVM)的詳細(xì)內(nèi)容，更多關(guān)于Matlab 圖像煙霧識(shí)別的資料請(qǐng)關(guān)注服務(wù)器之家其它相關(guān)文章！
原文鏈接：https://blog.csdn.net/qq_43215054/article/details/121734161
