人臉表情識別筆記（二）特徵提取之LBP（區域性二值模式）原理及MATLAB程式碼
一：原理部分
LBP（Local Binary Pattern，區域性二值模式）是一種用來描述影象區域性紋理特徵的運算元；它具有旋轉不變性和灰度不變性等顯著的優點。它是首先由T. Ojala, M.Pietikäinen, 和D. Harwood 在1994年提出，用於紋理特徵提取。而且，提取的特徵是影象的區域性的紋理特徵；
1、LBP特徵的描述
原始的LBP運算元定義為在3*3的視窗內，以視窗中心畫素為閾值，將相鄰的8個畫素的灰度值與其進行比較，若周圍畫素值大於中心畫素值，則該畫素點的位置被標記為1，否則為0。這樣，3*3鄰域內的8個點經比較可產生8位二進位制數（通常轉換為十進位制數即LBP
LBP的改進版本：
原始的LBP提出後，研究人員不斷對其提出了各種改進和優化。
（1）圓形LBP運算元：
基本的 LBP運算元的最大缺陷在於它只覆蓋了一個固定半徑範圍內的小區域，這顯然不能滿足不同尺寸和頻率紋理的需要。為了適應不同尺度的紋理特徵，並達到灰度和旋轉不變性的要求，Ojala等對 LBP 運算元進行了改進，將 3×3鄰域擴充套件到任意鄰域，並用圓形鄰域代替了正方形鄰域，改進後的 LBP 運算元允許在半徑為 R 的圓形鄰域內有任意多個畫素點。從而得到了諸如半徑為R
（2）LBP旋轉不變模式
從 LBP 的定義可以看出，LBP 運算元是灰度不變的，但卻不是旋轉不變的。影象的旋轉就會得到不同的 LBP值。
Maenpaa等人又將 LBP運算元進行了擴充套件，提出了具有旋轉不變性的 LBP 運算元，即不斷旋轉圓形鄰域得到一系列初始定義的 LBP值，取其最小值作為該鄰域的 LBP 值。
圖 2.5 給出了求取旋轉不變的 LBP 的過程示意圖，圖中運算元下方的數字表示該運算元對應的 LBP值，圖中所示的 8 種 LBP模式，經過旋轉不變的處理，最終得到的具有旋轉不變性的 LBP
（3）LBP等價模式
一個LBP運算元可以產生不同的二進位制模式，對於半徑為R的圓形區域內含有P個取樣點的LBP運算元將會產生P2種模式。很顯然，隨著鄰域集內取樣點數的增加，二進位制模式的種類是急劇增加的。例如：5×5鄰域內20個取樣點，有220＝1,048,576種二進位制模式。如此多的二值模式無論對於紋理的提取還是對於紋理的識別、分類及資訊的存取都是不利的。同時，過多的模式種類對於紋理的表達是不利的。例如，將LBP運算元用於紋理分類或人臉識別時，常採用LBP模式的統計直方圖來表達影象的資訊，而較多的模式種類將使得資料量過大，且直方圖過於稀疏。因此，需要對原始的LBP模式進行降維，使得資料量減少的情況下能最好的代表影象的資訊。
為了解決二進位制模式過多的問題，提高統計性，Ojala提出了採用一種“等價模式”（Uniform Pattern）來對LBP運算元的模式種類進行降維。Ojala等認為，在實際影象中，絕大多數LBP模式最多隻包含兩次從1到0或從0到1的跳變。因此，Ojala將“等價模式”定義為：當某個LBP所對應的迴圈二進位制數從0到1或從1到0最多有兩次跳變時，該LBP所對應的二進位制就稱為一個等價模式類。如00000000（0次跳變），00000111（只含一次從0到1的跳變），10001111（先由1跳到0，再由0跳到1，共兩次跳變）都是等價模式類。除等價模式類以外的模式都歸為另一類，稱為混合模式類，例如10010111（共四次跳變）（這是我的個人理解，不知道對不對）。
通過這樣的改進，二進位制模式的種類大大減少，而不會丟失任何資訊。模式數量由原來的2P種減少為 P ( P-1)+2種，其中P表示鄰域集內的取樣點數。對於3×3鄰域內8個取樣點來說，二進位制模式由原始的256種減少為58種，這使得特徵向量的維數更少，並且可以減少高頻噪聲帶來的影響。
2、LBP特徵用於檢測的原理
顯而易見的是，上述提取的LBP運算元在每個畫素點都可以得到一個LBP“編碼”，那麼，對一幅影象（記錄的是每個畫素點的灰度值）提取其原始的LBP運算元之後，得到的原始LBP特徵依然是“一幅圖片”（記錄的是每個畫素點的LBP值）。
LBP的應用中，如紋理分類、人臉分析等，一般都不將LBP圖譜作為特徵向量用於分類識別，而是採用LBP特徵譜的統計直方圖作為特徵向量用於分類識別。
因為，從上面的分析我們可以看出，這個“特徵”跟位置資訊是緊密相關的。直接對兩幅圖片提取這種“特徵”，並進行判別分析的話，會因為“位置沒有對準”而產生很大的誤差。後來，研究人員發現，可以將一幅圖片劃分為若干的子區域，對每個子區域內的每個畫素點都提取LBP特徵，然後，在每個子區域內建立LBP特徵的統計直方圖。如此一來，每個子區域，就可以用一個統計直方圖來進行描述；整個圖片就由若干個統計直方圖組成；
例如：一幅100*100畫素大小的圖片，劃分為10*10=100個子區域（可以通過多種方式來劃分區域），每個子區域的大小為10*10畫素；在每個子區域內的每個畫素點，提取其LBP特徵，然後，建立統計直方圖；這樣，這幅圖片就有10*10個子區域，也就有了10*10個統計直方圖，利用這10*10個統計直方圖，就可以描述這幅圖片了。之後，我們利用各種相似性度量函式，就可以判斷兩幅影象之間的相似性了；
3、對LBP特徵向量進行提取的步驟
（1）首先將檢測視窗劃分為16×16的小區域（cell）；
（2）對於每個cell中的一個畫素，將相鄰的8個畫素的灰度值與其進行比較，若周圍畫素值大於中心畫素值，則該畫素點的位置被標記為1，否則為0。這樣，3*3鄰域內的8個點經比較可產生8位二進位制數，即得到該視窗中心畫素點的LBP值；
（3）然後計算每個cell的直方圖，即每個數字（假定是十進位制數LBP值）出現的頻率；然後對該直方圖進行歸一化處理。
（4）最後將得到的每個cell的統計直方圖進行連線成為一個特徵向量，也就是整幅圖的LBP紋理特徵向量；
然後便可利用SVM或者其他機器學習演算法進行分類了。
二、程式碼部分
說明：一共有三個m檔案，一個是lbp.m, 存放主要的lbp演算法，一個是getmapping,用以做演算法的輔助函式，一個是lbptest.m,存放著測試程式碼。這三個檔案需要放到同一個資料夾，並在資料夾中新增相應的圖片，具體的圖片名字見lbptest.m的程式碼，執行lbptest.m可以檢視結果。程式碼最後給出效果圖
lbp.m %LBP returns the local binary pattern image or LBP histogram of an image.% J = LBP(I,R,N,MAPPING,MODE) returns either a local binary pattern
% coded image or the local binary pattern histogram of an intensity
% image I. The LBP codes are computed using N sampling points on a
% circle of radius R and using mapping table defined by MAPPING.
% See the getmapping function for different mappings and use 0 for
% no mapping. Possible values for MODE are
% 'h' or 'hist' to get a histogram of LBP codes
% 'nh' to get a normalized histogram
% Otherwise an LBP code image is returned.
%
% J = LBP(I) returns the original (basic) LBP histogram of image I
%
% J = LBP(I,SP,MAPPING,MODE) computes the LBP codes using n sampling
% points defined in (n * 2) matrix SP. The sampling points should be
% defined around the origin (coordinates (0,0)).
%
% Examples
% --------
% I=imread('rice.png');
% mapping=getmapping(8,'u2');
% H1=LBP(I,1,8,mapping,'h'); %LBP histogram in (8,1) neighborhood
% %using uniform patterns
% subplot(2,1,1),stem(H1);
%
% H2=LBP(I);
% subplot(2,1,2),stem(H2);
%
% SP=[-1 -1; -1 0; -1 1; 0 -1; -0 1; 1 -1; 1 0; 1 1];
% I2=LBP(I,SP,0,'i'); %LBP code image using sampling points in SP
% %and no mapping. Now H2 is equal to histogram
% %of I2.
function result = lbp(varargin) % image,radius,neighbors,mapping,mode)
% Version 0.3.2
% Authors: Marko Heikkil?and Timo Ahonen
% Changelog
% Version 0.3.2: A bug fix to enable using mappings together with a
% predefined spoints array
% Version 0.3.1: Changed MAPPING input to be a struct containing the mapping
% table and the number of bins to make the function run faster with high number
% of sampling points. Lauge Sorensen is acknowledged for spotting this problem.
% 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
getmapping.m
function mapping = getmapping(samples,mappingtype)
% Version 0.1.1
% Authors: Marko Heikkil?and Timo Ahonen
% Changelog
% 0.1.1 Changed output to be a structure
% Fixed a bug causing out of memory errors when generating rotation
% invariant mappings with high number of sampling points.
% Lauge Sorensen is acknowledged for spotting this problem.
table = 0:2^samples-1;
newMax = 0; %number of patterns in the resulting LBP code
index = 0;
if strcmp(mappingtype,'u2') %Uniform 2
newMax = samples*(samples-1) + 3;
for i = 0:2^samples-1
j = bitset(bitshift(i,1,samples),1,bitget(i,samples)); %rotate left
numt = sum(bitget(bitxor(i,j),1:samples)); %number of 1->0 and
%0->1 transitions
%in binary string
%x is equal to the
%number of 1-bits in
%XOR(x,Rotate left(x))
if numt <= 2
table(i+1) = index;
index = index + 1;
else
table(i+1) = newMax - 1;
end
end
end
if strcmp(mappingtype,'ri') %Rotation invariant
tmpMap = zeros(2^samples,1) - 1;
for i = 0:2^samples-1
rm = i;
r = i;
for j = 1:samples-1
r = bitset(bitshift(r,1,samples),1,bitget(r,samples)); %rotate
%left
if r < rm
rm = r;
end
end
if tmpMap(rm+1) < 0
tmpMap(rm+1) = newMax;
newMax = newMax + 1;
end
table(i+1) = tmpMap(rm+1);
end
end
if strcmp(mappingtype,'riu2') %Uniform & Rotation invariant
newMax = samples + 2;
for i = 0:2^samples - 1
j = bitset(bitshift(i,1,samples),1,bitget(i,samples)); %rotate left
numt = sum(bitget(bitxor(i,j),1:samples));
if numt <= 2
table(i+1) = sum(bitget(i,1:samples));
else
table(i+1) = samples+1;
end
end
end
mapping.table=table;
mapping.samples=samples;
mapping.num=newMax;
ibptest.m
%%%%%%%%%%%%%%%%%%%LBP變換後的影象%%%%%%%%%%%%%%%%%%
I11 = imread('test1.bmp');
SP=[-1 -1; -1 0; -1 1; 0 -1; -0 1; 1 -1; 1 0; 1 1];
I12=LBP(I11,SP,0,'i');
imshow(I12);
I21 = imread('test2.bmp');
SP=[-1 -1; -1 0; -1 1; 0 -1; -0 1; 1 -1; 1 0; 1 1];
I22=LBP(I21,SP,0,'i');
re1 = abs(I22-I12);
% re1 = 255 - re1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
I31 = imread('test3.bmp');
SP=[-1 -1; -1 0; -1 1; 0 -1; -0 1; 1 -1; 1 0; 1 1];
I32=LBP(I31,SP,0,'i');
I41 = imread('test4.bmp');
SP=[-1 -1; -1 0; -1 1; 0 -1; -0 1; 1 -1; 1 0; 1 1];
I42=LBP(I41,SP,0,'i');
re2 = abs(I32-I42);
% re2 = 255 - re2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
re3 = abs(I32 - I12);
% re3 = 255 - abs(I32 - I12);
re4 = abs(I32 - I22);
re5 = abs(I42 - I12);
re6 = abs(I42 - I22);
%%%%%%%%%%%%%%%%%%%%%%直方圖%%%%%%%%%%%%%%%%%%%%%%%%
mapping=getmapping(8,'u2');
H11=LBP(I11,1,8,mapping,'h'); %LBP histogram in (8,1) neighborhood using uniform patterns
subplot(2,2,1),stem(H11);
H12=LBP(I11);
subplot(2,2,2),stem(H12);
H21=LBP(I21,1,8,mapping,'h'); %LBP histogram in (8,1) neighborhood using uniform patterns
subplot(2,2,3),stem(H21);
H22=LBP(I21);
subplot(2,2,4),stem(H22);
H31=LBP(I31,1,8,mapping,'h'); %LBP histogram in (8,1) neighborhood using uniform patterns
subplot(2,2,3),stem(H31);
H32=LBP(I31);
subplot(2,2,4),stem(H32);
%%%%%%%%%%%%%%%%%%%%%%%無聊瞎寫%%%%%%%%%%%%%%%%%%%%%%%%
[row, col] = size(I11);
It1 = ones(row, col);
for i = 2 : col
It1(:, i - 1) = abs(I11(:, i) - I11(:, i - 1));
end
It2 = ones(row, col);
for i = 2 : row
It2(i-1,:) = abs(I11(i,:) - I11(i-1,:));
end