如何利用Matlab仿真實現(xiàn)圖像煙霧識別

這篇文章主要為大家展示了“如何利用Matlab仿真實現(xiàn)圖像煙霧識別”，內(nèi)容簡而易懂，條理清晰，希望能夠幫助大家解決疑惑，下面讓小編帶領大家一起研究并學習一下“如何利用Matlab仿真實現(xiàn)圖像煙霧識別”這篇文章吧。

一、算法簡介

1.1 c-means聚類算法

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

也就是說，聚類的目標是得到較高的類內(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個像素的灰度值與其進行比較，若周圍像素值大于中心像素值，則該像素點的位置被標記為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的的學習策略就是間隔最大化，可形式化為一個求解凸二次規(guī)劃的問題，也等價于正則化的合頁損失函數(shù)的最小化問題。SVM的的學習算法就是求解凸二次規(guī)劃的最優(yōu)化算法。

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

圖1-1SVM算法示意圖

二、算法實現(xiàn)

2.1 煙霧識別算法流程

1）首先對所有圖像進行預處理，假定將有煙當作正樣本，將沒煙看作負樣本，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聚類算法對訓練集和測試集圖像的像素進行聚類，實現(xiàn)圖像分割。

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

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

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

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

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

圖2-3 算法流程框圖

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

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

本文對預處理過后的訓練集和測試集圖像進行像素聚類，在此分別列舉一張有煙圖和無煙圖的圖像分割前后的效果對比。如圖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算法提取的特征進行特征降維，在此取前兩維特征進行模型訓練，前兩維度保留的信息含有98.75%，如下圖2-8所示.

2.4 SVM算法實現(xiàn)

在經(jīng)過上述圖像預處理、圖像像素聚類、LBP特征提取、PCA特征降維至兩維過程之后，將二維特征向量作為輸入訓練SVM模型，最終得到模型在訓練集上的分類準確度。

利用k-means+LBP+PCA+SVM算法，多次訓練模型，最終取平均值，得到在訓練集上的分類準確度為79%，在測試集上的分類準確度為78%。下圖為模型在訓練集上的分類效果圖。

三、結(jié)果分析

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

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

這是兩種不一樣的解決問題的思路。若采用本文的思路是Pipeline，若直接采用LBP+SVM的思路叫做end2end，各有優(yōu)缺點。Pipeline是將一個問題拆解成若干個子問題一次解決，然后串在一起，這種方法易于實現(xiàn)，且靈活性和可解釋性更高，但缺點是多個子任務會造成錯誤累積。end2end是將一個問題看成一個整體，一般可以獲得比pipeline更高的性能，但是整體像一個黑盒，可解釋性差。現(xiàn)在深度學習最新研究的趨勢是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程序\模式識別大作業(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
%讀取訓練集負樣本并計算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
%讀取測試集負樣本并計算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(:,:));%訓練集數(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);   %隨機選擇訓練集合測試集
cp = classperf(label_train);  %評估分類器性能
svmStruct3hog = svmtrain(pcaData1(train,1:k),label_train(train));%訓練SVM分類器  
%使用svmtrain進行訓練,得到訓練后的結(jié)構(gòu)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的值即為分類結(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);%計算預測的正確率

%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)圖像煙霧識別”這篇文章的所有內(nèi)容，感謝各位的閱讀！相信大家都有了一定的了解，希望分享的內(nèi)容對大家有所幫助，如果還想學習更多知識，歡迎關注億速云行業(yè)資訊頻道！