一、算法簡介

1.1 c-means聚類算法

聚類分析是根據(jù)在數(shù)據(jù)中發(fā)現(xiàn)的描述對象及其關(guān)系的信息，將數(shù)據(jù)對象進行分組。目的是使組內(nèi)的對象相互之間是相似的（相關(guān)的），而不同組中的對象是不同的（不相關(guān)的）。組內(nèi)相似性越大，組間差距越大，說明聚類效果越好。

也就是說，聚類的目標(biāo)是得到較高的類內(nèi)相似度和較低的類間相似度，使得類間的距離盡可能大，類內(nèi)樣本與類中心的距離盡可能小。在此，我們選用k-means聚類算法。

1 .2 LBP算法

LBP（Local Binary Pattern，局部二值模式）是一種用來描述圖像局部紋理特征的算子；它具有旋轉(zhuǎn)不變性和灰度不變性等顯著的優(yōu)點。它是首先由T. Ojala, M.Pietikäinen, 和D. Harwood 在1994年提出，用于紋理特征提取，提取的特征是圖像的局部的紋理特征。

原始的LBP算子定義為在3*3的窗口內(nèi)，以窗口中心像素為閾值，將相鄰的8個像素的灰度值與其進行比較，若周圍像素值大于中心像素值，則該像素點的位置被標(biāo)記為1，否則為0。這樣，3*3鄰域內(nèi)的8個點經(jīng)比較可產(chǎn)生8位二進制數(shù)（通常轉(zhuǎn)換為十進制數(shù)即LBP碼，共256種），即得到該窗口中心像素點的LBP值，并用這個值來反映該區(qū)域的紋理信息。

1.3 PCA算法

PCA(Principal Component Analysis)，即主成分分析方法，是一種使用最廣泛的數(shù)據(jù)降維算法。其算法步驟如下：    

1）數(shù)據(jù)中心化——去均值，根據(jù)需要，有的需要歸一化——Normalized；

2）求解協(xié)方差矩陣；

3）利用特征值分解/奇異值分解 求解特征值以及特征向量；

4）將特征值從大到小排序，保留前k個特征向量

5）利用特征向量構(gòu)造投影矩陣；

6）利用投影矩陣，得出降維的數(shù)據(jù)。

1.4 SVM算法

支持向量機（support vector machines, SVM）是一種二分類模型，它的基本模型是定義在特征空間上的間隔最大的線性分類器，間隔最大使它有別于感知機；SVM還包括核技巧，這使它成為實質(zhì)上的非線性分類器。SVM的的學(xué)習(xí)策略就是間隔最大化，可形式化為一個求解凸二次規(guī)劃的問題，也等價于正則化的合頁損失函數(shù)的最小化問題。SVM的的學(xué)習(xí)算法就是求解凸二次規(guī)劃的最優(yōu)化算法。

SVM學(xué)習(xí)的基本想法是求解能夠正確劃分訓(xùn)練數(shù)據(jù)集并且?guī)缀伍g隔最大的分離超平面。如下圖所示即為分類超平面，對于線性可分的數(shù)據(jù)集來說，這樣的超平面有無窮多個（即感知機），但是幾何間隔最大的分類超平面卻是唯一的。如下圖1-1SVM算法示意圖

圖1-1SVM算法示意圖

二、算法實現(xiàn)

2.1 煙霧識別算法流程

1）首先對所有圖像進行預(yù)處理，假定將有煙當(dāng)作正樣本，將沒煙看作負(fù)樣本，train集的smoke文件夾改名為pos，train集的non文件夾改名為neg；同理將test集的smoke文件夾改名為pos，test集的non文件夾改名為neg。為了對所有圖片進行處理，將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聚類算法對訓(xùn)練集和測試集圖像的像素進行聚類，實現(xiàn)圖像分割。

3）利用LBP對分割后的訓(xùn)練集圖像和測試集圖像進行特征提取。

4）分別對訓(xùn)練集和測試集使用主成分分析法（PCA）進行特征降維。

5）利用對訓(xùn)練集降維后得到的二維特征訓(xùn)練SVM二分類模型，

6）最后利用對測試集降維后得到的二維特征進行分類預(yù)測。

整體算法流程如下圖2-3所示

圖2-3 算法流程框圖

2.2 c-means算法實現(xiàn)

圖像分割是利用圖像的灰度、顏色、紋理、形狀等特征，把圖像分成若干個互不重疊的區(qū)域，并使這些特征在同一區(qū)域內(nèi)呈現(xiàn)相似性，在不同的區(qū)域之間存在明顯的差異性。然后就可以將分割的圖像中具有獨特性質(zhì)的區(qū)域提取出來用于不同的研究。圖像識別的基礎(chǔ)是圖像分割，其作用是把反映物體真實情況的、占據(jù)不同區(qū)域的、具有不同特性的目標(biāo)區(qū)分開來，并形成數(shù)字特征。因此本文利用c-means聚類算法實現(xiàn)圖像分割，實現(xiàn)對噪聲的過濾，在構(gòu)建煙霧識別模型的過程中，首先分別對無煙和有煙的圖像進行c-means聚類圖像分割。

本文對預(yù)處理過后的訓(xùn)練集和測試集圖像進行像素聚類，在此分別列舉一張有煙圖和無煙圖的圖像分割前后的效果對比。如圖2-4和圖2-5所示

圖2-4 無煙圖像分割前后對照圖

圖2-5有煙圖像分割前后對照圖 

2.3 LBP算法實現(xiàn)

本文LBP算法將像素聚類（3類）以后的圖像進行特征提取。在此分別列舉一張有煙圖和無煙圖的圖像特征提取前后的效果對比。

圖2-6無煙三像素聚類LBP特征提取前后對照圖 

圖2-7有煙三像素聚類LBP特征提取前后對照圖

本文PCA算法將HOG或LBP提取的特征進行特征降維，使數(shù)據(jù)可視化。PCA算法可以獲取原有特征的大部分信息，降維以后的前k個特征值保留下來的信息占原有信息的比例可有下式計算獲得。

對LBP算法提取的特征進行特征降維，在此取前兩維特征進行模型訓(xùn)練，前兩維度保留的信息含有98.75%，如下圖2-8所示.

2.4 SVM算法實現(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%，在測試集上的分類準(zhǔn)確度為78%。下圖為模型在訓(xùn)練集上的分類效果圖。

三、結(jié)果分析

經(jīng)過第二章的算法實現(xiàn)，最終得到了完整的SVM二分類模型，利用該模型對test中的pos樣本的圖片和neg樣本的圖片進行預(yù)測。預(yù)測前，首先需要對測試集圖片經(jīng)過預(yù)處理、其次利用k-means3聚類法對像素進行聚類得到最終圖像分割聚類圖、然后對聚類圖進行LBP特征提取、最后再利用PCA對提取出來的特征進行特征降維。將最終得到的二維特征向量作為模型的輸入，進行分類預(yù)測，最終得到結(jié)果。對于LBP特征提取方法，在訓(xùn)練集和測試集上的準(zhǔn)確率分別為79%和78%。經(jīng)過對比可以發(fā)現(xiàn)模型的泛化性能良好。

最后筆者不得不提的是，之所以采取上訴方法實現(xiàn)煙霧識別是因為，大作業(yè)要求必須包含聚類、分類、降維。筆者也嘗試過直接使用LBP+SVM實現(xiàn)煙霧識別的方法，并且對測試集的準(zhǔn)確率可以達到93%。

這是兩種不一樣的解決問題的思路。若采用本文的思路是Pipeline，若直接采用LBP+SVM的思路叫做end2end，各有優(yōu)缺點。Pipeline是將一個問題拆解成若干個子問題一次解決，然后串在一起，這種方法易于實現(xiàn)，且靈活性和可解釋性更高，但缺點是多個子任務(wù)會造成錯誤累積。end2end是將一個問題看成一個整體，一般可以獲得比pipeline更高的性能，但是整體像一個黑盒，可解釋性差。現(xiàn)在深度學(xué)習(xí)最新研究的趨勢是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');%載入測試正樣本列表  
sz_pos=size(ReadList_pos);   
label_pos=ones(sz_pos(1),1); %正樣本標(biāo)簽
ReadList_neg  = textread('neg_test_list.txt','%s','delimiter','\n');%載入測試負(fù)樣本列表
sz_neg=size(ReadList_neg);  
label_neg=zeros(sz_neg(1),1);%負(fù)樣本標(biāo)簽  
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);
%% 提取特征
%讀取訓(xùn)練集正樣本并計算lbp特征 
for i=1:sz1(1)
 name=char(ReadList1(i,1));  
 image1=imread(strcat('F:\模式識別matlab程序\模式識別大作業(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ù)樣本并計算lbp特征  
for j=1:sz2(1)
 name= char(ReadList2(j,1));  
 image2=imread(strcat('F:\模式識別matlab程序\模式識別大作業(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
%讀取測試集正樣本并計算lbp特征
for m=1:sz_pos(1)
 test_name= char(ReadList_pos(m,1));  
 image3=imread(strcat('F:\模式識別matlab程序\模式識別大作業(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
%讀取測試集負(fù)樣本并計算lbp特征
for n =1:sz_neg(1)
  test_name=char(ReadList_neg(n,1)); 
  image4=imread(strcat('F:\模式識別matlab程序\模式識別大作業(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個特征值所包含的圖像信息比例
x0 = bsxfun(@minus,data2,mean(data2,1));
pcaData2_sw = x0*COEFF(:,:);
pcaData2 = pcaData2_sw(:,1:k);
%%  評估方法：交叉驗證法
[train, test] = crossvalind('holdOut',label_train);   %隨機選擇訓(xùn)練集合測試集
cp = classperf(label_train);  %評估分類器性能
svmStruct3hog = svmtrain(pcaData1(train,1:k),label_train(train));%訓(xùn)練SVM分類器  
%使用svmtrain進行訓(xùn)練,得到訓(xùn)練后的結(jié)構(gòu)svmStruct3hog,在預(yù)測時使用
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的值即為分類結(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);%計算預(yù)測的正確率

%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ù)表示隨機選擇K個聚類中心
% F:樣本聚類編號
if nargin~=2
  error('IMKMEANS:InputParamterNotRight','只能有兩個輸入?yún)?shù)！');
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.循環(huán)搜索聚類中心
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仿真實現(xiàn)圖像煙霧識別(k-means聚類圖像分割+LBP+PCA+SVM)的詳細(xì)內(nèi)容，更多關(guān)于Matlab 圖像煙霧識別的資料請關(guān)注服務(wù)器之家其它相關(guān)文章！

原文鏈接：https://blog.csdn.net/qq_43215054/article/details/121734161