一、演算法簡介
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更高的性能，但是整體像一個黑盒，可解釋性差，現在深度學習最新研究的趨勢是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]);
轉載請註明出處，本文鏈接：https://www.uj5u.com/qita/374542.html
標籤：AI