一、算法簡介
1.1 c-means聚類算法
聚類分析是根據在數據中發現的描述對象及其關系的信息，將數據對象進行分組。目的是使組內的對象相互之間是相似的（相關的），而不同組中的對象是不同的（不相關的）。組內相似性越大，組間差距越大，說明聚類效果越好。
也就是說，聚類的目標是得到較高的類內相似度和較低的類間相似度，使得類間的距離盡可能大，類內樣本與類中心的距離盡可能小。在此，我們選用k-means聚類算法。
1 .2 LBP算法
LBP（Local Binary Pattern，局部二值模式）是一種用來描述圖像局部紋理特征的算子；它具有旋轉不變性和灰度不變性等顯著的優點。它是首先由T. Ojala, M.Pietikäinen, 和D. Harwood 在1994年提出，用于紋理特征提取，提取的特征是圖像的局部的紋理特征。
原始的LBP算子定義為在3*3的窗口內，以窗口中心像素為閾值，將相鄰的8個像素的灰度值與其進行比較，若周圍像素值大于中心像素值，則該像素點的位置被標記為1，否則為0。這樣，3*3鄰域內的8個點經比較可產生8位二進制數（通常轉換為十進制數即LBP碼，共256種），即得到該窗口中心像素點的LBP值，并用這個值來反映該區域的紋理信息。
1.3 PCA算法
PCA(Principal Component Analysis)，即主成分分析方法，是一種使用最廣泛的數據降維算法。其算法步驟如下：
1）數據中心化——去均值，根據需要，有的需要歸一化——Normalized；
2）求解協方差矩陣；
3）利用特征值分解/奇異值分解 求解特征值以及特征向量；
4）將特征值從大到小排序，保留前k個特征向量
5）利用特征向量構造投影矩陣；
6）利用投影矩陣，得出降維的數據。
1.4 SVM算法
支持向量機（support vector machines, SVM）是一種二分類模型，它的基本模型是定義在特征空間上的間隔最大的線性分類器，間隔最大使它有別于感知機；SVM還包括核技巧，這使它成為實質上的非線性分類器。SVM的的學習策略就是間隔最大化，可形式化為一個求解凸二次規劃的問題，也等價于正則化的合頁損失函數的最小化問題。SVM的的學習算法就是求解凸二次規劃的最優化算法。
SVM學習的基本想法是求解能夠正確劃分訓練數據集并且幾何間隔最大的分離超平面。如下圖所示即為分類超平面，對于線性可分的數據集來說，這樣的超平面有無窮多個（即感知機），但是幾何間隔最大的分類超平面卻是唯一的。如下圖1-1SVM算法示意圖
圖1-1SVM算法示意圖
二、算法實現
2.1 煙霧識別算法流程
1）首先對所有圖像進行預處理，假定將有煙當作正樣本，將沒煙看作負樣本，train集的smoke文件夾改名為pos，train集的non文件夾改名為neg；同理將test集的smoke文件夾改名為pos，test集的non文件夾改名為neg。為了對所有圖片進行處理，將train和test中的pos和neg中的圖片全部規范命名格式為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聚類算法對訓練集和測試集圖像的像素進行聚類，實現圖像分割。
3）利用LBP對分割后的訓練集圖像和測試集圖像進行特征提取。
4）分別對訓練集和測試集使用主成分分析法（PCA）進行特征降維。
5）利用對訓練集降維后得到的二維特征訓練SVM二分類模型，
6）最后利用對測試集降維后得到的二維特征進行分類預測。
整體算法流程如下圖2-3所示
圖2-3 算法流程框圖
2.2 c-means算法實現
圖像分割是利用圖像的灰度、顏色、紋理、形狀等特征，把圖像分成若干個互不重疊的區域，并使這些特征在同一區域內呈現相似性，在不同的區域之間存在明顯的差異性。然后就可以將分割的圖像中具有獨特性質的區域提取出來用于不同的研究。圖像識別的基礎是圖像分割，其作用是把反映物體真實情況的、占據不同區域的、具有不同特性的目標區分開來，并形成數字特征。因此本文利用c-means聚類算法實現圖像分割，實現對噪聲的過濾，在構建煙霧識別模型的過程中，首先分別對無煙和有煙的圖像進行c-means聚類圖像分割。
本文對預處理過后的訓練集和測試集圖像進行像素聚類，在此分別列舉一張有煙圖和無煙圖的圖像分割前后的效果對比。如圖2-4和圖2-5所示
圖2-4 無煙圖像分割前后對照圖
圖2-5有煙圖像分割前后對照圖
2.3 LBP算法實現
本文LBP算法將像素聚類（3類）以后的圖像進行特征提取。在此分別列舉一張有煙圖和無煙圖的圖像特征提取前后的效果對比。
圖2-6無煙三像素聚類LBP特征提取前后對照圖
圖2-7有煙三像素聚類LBP特征提取前后對照圖
本文PCA算法將HOG或LBP提取的特征進行特征降維，使數據可視化。PCA算法可以獲取原有特征的大部分信息，降維以后的前k個特征值保留下來的信息占原有信息的比例可有下式計算獲得。
對LBP算法提取的特征進行特征降維，在此取前兩維特征進行模型訓練，前兩維度保留的信息含有98.75%，如下圖2-8所示.
2.4 SVM算法實現
在經過上述圖像預處理、圖像像素聚類、LBP特征提取、PCA特征降維至兩維過程之后，將二維特征向量作為輸入訓練SVM模型，最終得到模型在訓練集上的分類準確度。
利用k-means+LBP+PCA+SVM算法，多次訓練模型，最終取平均值，得到在訓練集上的分類準確度為79%，在測試集上的分類準確度為78%。下圖為模型在訓練集上的分類效果圖。
三、結果分析
經過第二章的算法實現，最終得到了完整的SVM二分類模型，利用該模型對test中的pos樣本的圖片和neg樣本的圖片進行預測。預測前，首先需要對測試集圖片經過預處理、其次利用k-means3聚類法對像素進行聚類得到最終圖像分割聚類圖、然后對聚類圖進行LBP特征提取、最后再利用PCA對提取出來的特征進行特征降維。將最終得到的二維特征向量作為模型的輸入，進行分類預測，最終得到結果。對于LBP特征提取方法，在訓練集和測試集上的準確率分別為79%和78%。經過對比可以發現模型的泛化性能良好。
最后筆者不得不提的是，之所以采取上訴方法實現煙霧識別是因為，大作業要求必須包含聚類、分類、降維。筆者也嘗試過直接使用LBP+SVM實現煙霧識別的方法，并且對測試集的準確率可以達到93%。
這是兩種不一樣的解決問題的思路。若采用本文的思路是Pipeline，若直接采用LBP+SVM的思路叫做end2end，各有優缺點。Pipeline是將一個問題拆解成若干個子問題一次解決，然后串在一起，這種方法易于實現，且靈活性和可解釋性更高，但缺點是多個子任務會造成錯誤累積。end2end是將一個問題看成一個整體，一般可以獲得比pipeline更高的性能，但是整體像一個黑盒，可解釋性差?，F在深度學習最新研究的趨勢是end2end的方法。
%基于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程序\模式識別大作業\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
%讀取訓練集負樣本并計算lbp特征
for j=1:sz2(1)
name= char(ReadList2(j,1));
image2=imread(strcat('F:\模式識別matlab程序\模式識別大作業\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
%讀取測試集正樣本并計算lbp特征
for m=1:sz_pos(1)
test_name= char(ReadList_pos(m,1));
image3=imread(strcat('F:\模式識別matlab程序\模式識別大作業\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
%讀取測試集負樣本并計算lbp特征
for n =1:sz_neg(1)
test_name=char(ReadList_neg(n,1));
image4=imread(strcat('F:\模式識別matlab程序\模式識別大作業\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
%數據降維
[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);%計算預測的正確率
%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:聚類中心，可以是整數或者數組，整數表示隨機選擇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]);
以上就是利用Matlab仿真實現圖像煙霧識別(k-means聚類圖像分割+LBP+PCA+SVM)的詳細內容，更多關于Matlab 圖像煙霧識別的資料請關注服務器之家其它相關文章！
原文鏈接：https://blog.csdn.net/qq_43215054/article/details/121734161
