LBP特征提取

与第一篇博文特征脸方法不同，LBP（Local Binary Patterns，局部二值模式）是提取局部特征作为判别依据的。LBP方法显著的优点是对光照不敏感，但是依然没有解决姿态和表情的问题。不过相比于特征脸方法，LBP的识别率已经有了很大的提升。在[1]的文章里，有些人脸库的识别率已经达到了98%+。

1、LBP特征提取

最初的LBP是定义在像素3x3邻域内的，以邻域中心像素为阈值，将相邻的8个像素的灰度值与其进行比较，若周围像素值大于中心像素值，则该像素点的位置被标记为1，否则为0。这样，3x3邻域内的8个点经比较可产生8位二进制数（通常转换为十进制数即LBP码，共256种），即得到该邻域中心像素点的LBP值，并用这个值来反映该区域的纹理信息。如下图所示：

用比较正式的公式来定义的话：

其中代表3x3邻域的中心元素，它的像素值为ic，ip代表邻域内其他像素的值。s(x)是符号函数，定义如下：

LBP的改进版本

（1）圆形LBP算子

基本的 LBP算子的最大缺陷在于它只覆盖了一个固定半径范围内的小区域，这显然不能满足不同尺寸和频率纹理的需要。为了适应不同尺度的纹理特征，并达到灰度和旋转不变性的要求，Ojala等对 LBP 算子进行了改进，将 3×3邻域扩展到任意邻域，并用圆形邻域代替了正方形邻域，改进后的 LBP 算子允许在半径为 R 的圆形邻域内有任意多个像素点。从而得到了诸如半径为R的圆形区域内含有P个采样点的LBP算子。比如下图定了一个5x5的邻域：

上图内有八个黑色的采样点，每个采样点的值可以通过下式计算：

其中为邻域中心点，为某个采样点。通过上式可以计算任意个采样点的坐标，但是计算得到的坐标未必完全是整数，所以可以通过双线性插值来得到该采样点的像素值：

（2）LBP等价模式

一个LBP算子可以产生不同的二进制模式，对于半径为R的圆形区域内含有P个采样点的LBP算子将会产生2^P种模式。很显然，随着邻域集内采样点数的增加，二进制模式的种类是急剧增加的。例如：5×5邻域内20个采样点，有2²⁰＝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（共四次跳变）。比如下图给出了几种等价模式的示意图。

       通过这样的改进，二进制模式的种类大大减少，而不会丢失任何信息。模式数量由原来的2^P种减少为 P ( P-1)+2种，其中P表示邻域集内的采样点数。对于3×3邻域内8个采样点来说，二进制模式由原始的256种减少为58种，这使得特征向量的维数更少，并且可以减少高频噪声带来的影响。这几段摘自[2]。

通过上述方法，每个像素都会根据邻域信息得到一个LBP值，如果以图像的形式显示出来可以得到下图，明显LBP对光照有较强的鲁棒性。

2、LBP特征匹配

如果将以上得到的LBP值直接用于人脸识别，其实和不提取LBP特征没什么区别，会造成计算量准确率等一系列问题。文献[1]中，将一副人脸图像分为7x7的子区域（如下图），并在子区域内根据LBP值统计其直方图，以直方图作为其判别特征。这样做的好处是在一定范围内避免图像没完全对准的情况，同时也对LBP特征做了降维处理。

对于得到的直方图特征，有多种方法可以判别其相似性，假设已知人脸直方图为M_i​，待匹配人脸直方图为S_i，那么可以通过:

(1)直方图交叉核方法

该方法的介绍在博文：Histogram intersection(直方图交叉核,Pyramid Match Kernel)

(2)卡方统计方法

该方法的介绍在博文：卡方检验（Chi square statistic）

参考文献：

[1]Timo Ahonen, Abdenour Hadid：Face Recognition with Local Binary Patterns

[2]目标检测的图像特征提取之（二）LBP特征

图像分析：LBP特征解析与代码

分类： 算法与理论研究 图像分析专栏 2013-08-29 11:00  542人阅读  评论(0)  收藏  举报

LBP

LBP(Local Binary Patterns)，即局部二值模式，是一种描述图像局部空间结构的非参数算子。芬兰Oulu大学的T.Ojala等人于1996年提出这个算子用来分析图像纹理特征，并且描述了它在纹理分类中的强区分能力。LBP算子定义为一种灰度尺度不变的纹理算子，是从局部邻域纹理的普通定义得来的。

基本思想是：用中心像素的灰度值作为阈值，与它的邻域相比较得到的二进制码来表述局部纹理特征。

在纹理分析方面，LBP算子是最好的纹理描述符之一，它的主要优点有以下几点：

Ø 通过它的定义可知，LBP算子的灰度尺度不随任何单一变换而变化，因此灰度尺度的鲁棒性好，也就是光照条件下的鲁棒性好；

Ø 计算速度快。由于它可以通过在小邻域内进行比较操作得到，使得在复杂的实时条件下分析图像成为可能；

Ø 由于 LBP 算子是一种无参数 (Non-Parametric) 的方法，在应用过程中不需要对它的分布进行预先假设。

_____________________________________________________________________________________

[plain]  view plain copy

1. %LBP returns the local binary pattern image or LBP histogram of an image.  
2. %  J = LBP(I,R,N,MAPPING,MODE) returns either a local binary pattern  
3. %  coded image or the local binary pattern histogram of an intensity  
4. %  image I. The LBP codes are computed using N sampling points on a   
5. %  circle of radius R and using mapping table defined by MAPPING.   
6. %  See the getmapping function for different mappings and use 0 for  
7. %  no mapping. Possible values for MODE are  
8. %       'h' or 'hist'  to get a histogram of LBP codes  
9. %       'nh'           to get a normalized histogram  
10. %  Otherwise an LBP code image is returned.  
11. %  
12. %  J = LBP(I) returns the original (basic) LBP histogram of image I  
13. %  
14. %  J = LBP(I,SP,MAPPING,MODE) computes the LBP codes using n sampling  
15. %  points defined in (n * 2) matrix SP. The sampling points should be  
16. %  defined around the origin (coordinates (0,0)).  
17. %  
18. %  Examples  
19. %  --------  
20. %       I=imread('rice.png');  
21. %       mapping=getmapping(8,'u2');   
22. %       H1=LBP(I,1,8,mapping,'h'); %LBP histogram in (8,1) neighborhood  
23. %                                  %using uniform patterns  
24. %       subplot(2,1,1),stem(H1);  
25. %  
26. %       H2=LBP(I);  
27. %       subplot(2,1,2),stem(H2);  
28. %  
29. %       SP=[-1 -1; -1 0; -1 1; 0 -1; -0 1; 1 -1; 1 0; 1 1];  
30. %       I2=LBP(I,SP,0,'i'); %LBP code image using sampling points in SP  
31. %                           %and no mapping. Now H2 is equal to histogram  
32. %                           %of I2.  
33.   
34. function result = lbp(varargin) % image,radius,neighbors,mapping,mode)  
35. % Version 0.3.2  
36. % Authors: Marko Heikkil�and Timo Ahonen  
37.   
38. % Changelog  
39. % Version 0.3.2: A bug fix to enable using mappings together with a  
40. % predefined spoints array  
41. % Version 0.3.1: Changed MAPPING input to be a struct containing the mapping  
42. % table and the number of bins to make the function run faster with high number  
43. % of sampling points. Lauge Sorensen is acknowledged for spotting this problem.  
44.   
45.   
46. % Check number of input arguments.  
47. error(nargchk(1,5,nargin));  
48.   
49. image=varargin{1};  
50. d_image=double(image);  
51.   
52. if nargin==1  
53.     spoints=[-1 -1; -1 0; -1 1; 0 -1; -0 1; 1 -1; 1 0; 1 1];  
54.     neighbors=8;  
55.     mapping=0;  
56.     mode='h';  
57. end  
58.   
59. if (nargin == 2) && (length(varargin{2}) == 1)  
60.     error('Input arguments');  
61. end  
62.   
63. if (nargin > 2) && (length(varargin{2}) == 1)  
64.     radius=varargin{2};  
65.     neighbors=varargin{3};  
66.       
67.     spoints=zeros(neighbors,2);  
68.   
69.     % Angle step.  
70.     a = 2*pi/neighbors;  
71.       
72.     for i = 1:neighbors  
73.         spoints(i,1) = -radius*sin((i-1)*a);  
74.         spoints(i,2) = radius*cos((i-1)*a);  
75.     end  
76.       
77.     if(nargin >= 4)  
78.         mapping=varargin{4};  
79.         if(isstruct(mapping) && mapping.samples ~= neighbors)  
80.             error('Incompatible mapping');  
81.         end  
82.     else  
83.         mapping=0;  
84.     end  
85.       
86.     if(nargin >= 5)  
87.         mode=varargin{5};  
88.     else  
89.         mode='h';  
90.     end  
91. end  
92.   
93. if (nargin > 1) && (length(varargin{2}) > 1)  
94.     spoints=varargin{2};  
95.     neighbors=size(spoints,1);  
96.       
97.     if(nargin >= 3)  
98.         mapping=varargin{3};  
99.         if(isstruct(mapping) && mapping.samples ~= neighbors)  
100.             error('Incompatible mapping');  
101.         end  
102.     else  
103.         mapping=0;  
104.     end  
105.       
106.     if(nargin >= 4)  
107.         mode=varargin{4};  
108.     else  
109.         mode='h';  
110.     end     
111. end  
112.   
113. % Determine the dimensions of the input image.  
114. [ysize xsize] = size(image);  
115.   
116.   
117.   
118. miny=min(spoints(:,1));  
119. maxy=max(spoints(:,1));  
120. minx=min(spoints(:,2));  
121. maxx=max(spoints(:,2));  
122.   
123. % Block size, each LBP code is computed within a block of size bsizey*bsizex  
124. bsizey=ceil(max(maxy,0))-floor(min(miny,0))+1;  
125. bsizex=ceil(max(maxx,0))-floor(min(minx,0))+1;  
126.   
127. % Coordinates of origin (0,0) in the block  
128. origy=1-floor(min(miny,0));  
129. origx=1-floor(min(minx,0));  
130.   
131. % Minimum allowed size for the input image depends  
132. % on the radius of the used LBP operator.  
133. if(xsize < bsizex || ysize < bsizey)  
134.   error('Too small input image. Should be at least (2*radius+1) x (2*radius+1)');  
135. end  
136.   
137. % Calculate dx and dy;  
138. dx = xsize - bsizex;  
139. dy = ysize - bsizey;  
140.   
141. % Fill the center pixel matrix C.  
142. C = image(origy:origy+dy,origx:origx+dx);  
143. d_C = double(C);  
144.   
145. bins = 2^neighbors;  
146.   
147. % Initialize the result matrix with zeros.  
148. result=zeros(dy+1,dx+1);  
149.   
150. %Compute the LBP code image  
151.   
152. for i = 1:neighbors  
153.   y = spoints(i,1)+origy;  
154.   x = spoints(i,2)+origx;  
155.   % Calculate floors, ceils and rounds for the x and y.  
156.   fy = floor(y); cy = ceil(y); ry = round(y);  
157.   fx = floor(x); cx = ceil(x); rx = round(x);  
158.   % Check if interpolation is needed.  
159.   if (abs(x - rx) < 1e-6) && (abs(y - ry) < 1e-6)  
160.     % Interpolation is not needed, use original datatypes  
161.     N = image(ry:ry+dy,rx:rx+dx);  
162.     D = N >= C;   
163.   else  
164.     % Interpolation needed, use double type images   
165.     ty = y - fy;  
166.     tx = x - fx;  
167.   
168.     % Calculate the interpolation weights.  
169.     w1 = (1 - tx) * (1 - ty);  
170.     w2 =      tx  * (1 - ty);  
171.     w3 = (1 - tx) *      ty ;  
172.     w4 =      tx  *      ty ;  
173.     % Compute interpolated pixel values  
174.     N = w1*d_image(fy:fy+dy,fx:fx+dx) + w2*d_image(fy:fy+dy,cx:cx+dx) + ...  
175.         w3*d_image(cy:cy+dy,fx:fx+dx) + w4*d_image(cy:cy+dy,cx:cx+dx);  
176.     D = N >= d_C;   
177.   end    
178.   % Update the result matrix.  
179.   v = 2^(i-1);  
180.   result = result + v*D;  
181. end  
182.   
183. %Apply mapping if it is defined  
184. if isstruct(mapping)  
185.     bins = mapping.num;  
186.     for i = 1:size(result,1)  
187.         for j = 1:size(result,2)  
188.             result(i,j) = mapping.table(result(i,j)+1);  
189.         end  
190.     end  
191. end  
192.   
193. if (strcmp(mode,'h') || strcmp(mode,'hist') || strcmp(mode,'nh'))  
194.     % Return with LBP histogram if mode equals 'hist'.  
195.     result=hist(result(:),0:(bins-1));  
196.     if (strcmp(mode,'nh'))  
197.         result=result/sum(result);  
198.     end  
199. else  
200.     %Otherwise return a matrix of unsigned integers  
201.     if ((bins-1)<=intmax('uint8'))  
202.         result=uint8(result);  
203.     elseif ((bins-1)<=intmax('uint16'))  
204.         result=uint16(result);  
205.     else  
206.         result=uint32(result);  
207.     end  
208. end  
209.   
210. end  

[plain]  view plain copy

1. %GETMAPPING returns a structure containing a mapping table for LBP codes.  
2. %  MAPPING = GETMAPPING(SAMPLES,MAPPINGTYPE) returns a structure containing a mapping table for  
3. %  LBP codes in a neighbourhood of SAMPLES sampling  points. Possible values for MAPPINGTYPE are  
4. %       'u2'   for uniform LBP  
5. %       'ri'   for rotation-invariant LBP  
6. %       'riu2' for uniform rotation-invariant LBP.  
7. %  
8. %  Example:  
9. %       I=imread('rice.tif');  
10. %       MAPPING=getmapping(16,'riu2');  
11. %       LBPHIST=lbp(I,2,16,MAPPING,'hist');  
12. %  Now LBPHIST contains a rotation-invariant uniform LBP  histogram in a (16,2) neighbourhood.  
13. %  
14.   
15. function mapping = getmapping(samples,mappingtype)  
16. % Version 1.0 beta  
17. % Authors: Marko Heikkil?and Timo Ahonen  
18. % Editor： LinJianmin - Huaqiao University  
19. % Email:    jianmin1990@outlook.com  
20.   
21. table = 0:2^samples-1;  
22. newMax  = 0; %number of patterns in the resulting LBP code  
23. index   = 0;  
24.   
25. if strcmp(mappingtype,'u2') %Uniform 2  
26.   newMax = samples*(samples-1) + 3;   
27.     
28.   for i = 0:2^samples-1  
29.     j = bitset(bitshift(i,1,samples),1,bitget(i,samples)); %rotate left  
30.     numt = sum(bitget(bitxor(i,j),1:samples)); %number of 1->0 and 0->1 transitions in binary string   
31.                                                %x is equal to the number of 1-bits in XOR(x,Rotate left(x))   
32.     if numt <= 2  
33.       table(i+1) = index;  
34.       index = index + 1;  
35.     else  
36.       table(i+1) = newMax - 1;  
37.     end  
38.   end  
39. end  
40.   
41. if strcmp(mappingtype,'ri') %Rotation invariant  
42.   tmpMap = zeros(2^samples,1) - 1;  
43.   for i = 0:2^samples-1  
44.     rm = i;  
45.     r  = i;  
46.     for j = 1:samples-1  
47.       r = bitset(bitshift(r,1,samples),1,bitget(r,samples)); %rotate  
48.                                                              %left  
49.       if r < rm  
50.         rm = r;  
51.       end  
52.     end  
53.     if tmpMap(rm+1) < 0  
54.       tmpMap(rm+1) = newMax;  
55.       newMax = newMax + 1;  
56.     end  
57.     table(i+1) = tmpMap(rm+1);  
58.   end  
59. end  
60.   
61. if strcmp(mappingtype,'riu2') %Uniform & Rotation invariant  
62.   newMax = samples + 2;  
63.   for i = 0:2^samples - 1  
64.     j = bitset(bitshift(i,1,samples),1,bitget(i,samples)); %rotate left  
65.     numt = sum(bitget(bitxor(i,j),1:samples));  
66.     if numt <= 2  
67.       table(i+1) = sum(bitget(i,1:samples));  
68.     else  
69.       table(i+1) = samples+1;  
70.     end  
71.   end  
72. end  
73.   
74. mapping.table=table;  
75. mapping.samples=samples;  
76. mapping.num=newMax;