LBP特征

LBP 算子是一种有效的纹理描述算子, 它具有旋转不变性和灰度不变性等显著优点。其基本思想是用其中心像素的灰度值作为阈值,与它的邻域相比较得到的二进制码来表述局部纹理特征。

基本LBP

首先介绍LBP 算子的计算方式。对于一幅图像中的某个局部区域内的任意像素f(xc , yc ),以其为中心点gc , 对3 ×3 窗口内的8 个点g0 ,…, g7 , 纹理T 定义如下:

以窗口中心点灰度值为阈值对窗口内其它像素做二值化处理, 如公式(2)所示:

得到一个8 位的二进制数, 按公式(3)对像素不同位置进行加权求和, 即可得到该窗口的LBP 值:

一个基本的LBP 算子如下图所示。

基本LBP算子示意图

改进LBP

为了适应不同尺度的纹理特征,对LBP 算子进行了改进, 将3 ×3 邻域扩展到任意邻域, 并用圆形邻域代替了正方形邻域, 采用双线性插值算法计算没有完全落在像素位置的点的灰度值。此外, 改进后的LBP 算子允许在半径为R 的圆形邻域内有任意多个像素点, 如图3 所示。符号表示在半径为R 的圆形邻域内有P 个像素点。

几种LBP算子

原始LBP与改进LBP的差别：

（1）原始定义于方形邻域上,而修改过的定义于圆形邻域中,从而更易于定义具有旋转不变性的纹理描述算子

（2）修改后的,算子中,其邻域中对角线上的像素值是通过线性内插得到的

uniform LBP

注意到随着采样点数的增加, 二进制模式的种类会急剧增加。如3 ×3 邻域内8 个采样点, 对应有2^8 种二进制模式;5 ×5 邻域内20 个采样点, 对应有2^20 种二进制模式;7 ×7 的邻域内36个采样点, 对应的二进制模式种类多达2^36种。显然, 如此多的二值模式对于提取纹理是不利的。为了解决二进制模式过多的问题, 提高统计性,利用所谓的“等价模式类”(uniform patterns)对LBP 进行了改进。当某个局部二进制模式首尾相连环的形所对应的循环二进制数从0 到1 或从1 到0 最多有两次跳变时, 该局部二进制模式所对应的二进制就称为一个等价模式类, 如00000000 ,11111111 , 10001111 都是等价模式类, 除等价模式类以外的模式都归为另一类, 称为混合模式类。这使得基本LBP算子的二进制模式的种类由原来的2^P 种减少为P(P-1)+3 种，例如P=8时等价模式类58个，混合模式类1个，共59种。其中，P 表示邻域集内的采样点数。通过这样的改进, 局部二值模式的种类大大减少。等价模式类占总模式中的绝大多数, 利用这些等价模式类和混合模式类的直方图, 可提取更能代表图像本质特性的特征。

LBP用于人脸识别

人脸的分类识别是通过匹配LBP 直方图序列特征的相似度来实现分类，首先,面部区域被分成许多小区域,并在此基础上生成多个直方图。然后,将这些直方图连接,生成一个空间特征得以增强、有效表征面部图像的直方图。计算出分块完备LBP 直方图特征的不相似度，然后根据最近邻准则进行分类。目前常用的几种直方图间距离计算方法有：1) 直方图相交距离 ; 2) 对数似然统计法 ; 3) Chi 平方统计法。采用Chi 平方统计法进行直方图间的不相似性度量，两个样本间的Chi距离值越小，表明两个样本越相似。

在Chi平方相似性度量过程中,为了实现对图像的识别,需要对测试集和训练集中的每一幅图像的增强直方图向量进行两两Chi平方距离计算,得到对应的平方距离矩阵。设训练集特征向量矩阵为A，大小为m*h1,其中m为训练集所包含的图像数,h1为特征向量维数；测试集特征向量矩阵为B，大小为n*h2，其中n为测试集所包含的图像数,h2为特征向量维数；则计算所得的平方距离矩阵为m*n维，Cij表示训练集中的图像i与测试集中的图像j之间的Chi平方距离。通过计算测试样本的分块完备LBP 直方图特征和所有训练样本的分块完备LBP 直方图特征的Chi距离值，然后根据最近邻准则，将测试样本归为距离值最小的那一类。

LBP特征提取的matlab代码

%getLBPFea.m
function [histLBP, MatLBP] = getLBPFea(I)
% 计算分区图像 I 的LBP特征,(8,2),uniform
%
% 输入：I --- 分区图像
%
% 返回值: MatLBP --- LBP响应矩阵
%               histLBP --- 1维行向量，LBP直方图

% 获得分块图像I的大小
[m n] = size(I);
rad = 2;
if (m <= 2*rad) || (n <= 2*rad)
    error('I is too small to compute LBP feature!');
end

MatLBP = zeros(m-2*rad, n-2*rad);

% 读入 LBP 映射(像素灰度与直方图收集箱索引的映射)
load Mat/LBPMap.mat;

for ii = 1+rad : m-rad
    for jj = 1+rad : n-rad
        nCnt = 1;
        
        
        % 计算(8,2)邻域的像素值，不在像素中心的点通过双线性插值获得其值
        nbPT(nCnt) = I(ii, jj-rad);
        nCnt = nCnt + 1;
        
        horInterp1 = I(ii-2, jj-2) + 0.5858*( I(ii-2, jj-1) - I(ii-2, jj-2) ); % 水平方向插值
        horInterp2 = I(ii-1, jj-2) + 0.5858*( I(ii-1, jj-1) - I(ii-1, jj-2) ); % 水平方向插值
        verInterp = horInterp1 + 0.5858*( horInterp2 - horInterp1 ); % 竖直方向插值
        nbPT(nCnt) = verInterp;
        nCnt = nCnt + 1;
        
        nbPT(nCnt) = I(ii-2, jj);
        nCnt = nCnt + 1;
        
        horInterp1 = I(ii-2, jj+1) + 0.4142*( I(ii-2, jj+2) - I(ii-2, jj+1) );
        horInterp2 = I(ii-1, jj+1) + 0.4142*( I(ii-1, jj+2) - I(ii-1, jj+1) );
        verInterp = horInterp1 + 0.5858*( horInterp2 - horInterp1 );
        nbPT(nCnt) = verInterp;
        nCnt = nCnt + 1;
        
        nbPT(nCnt) = I(ii, jj+2);
        nCnt = nCnt + 1;
        
        horInterp1 = I(ii+1, jj+1) + 0.4142*( I(ii+1, jj+2) - I(ii+1, jj+1) );
        horInterp2 = I(ii+2, jj+1) + 0.4142*( I(ii+2, jj+2) - I(ii+2, jj+1) );
        verInterp = horInterp1 + 0.4142*( horInterp2 - horInterp1 );
        nbPT(nCnt) = verInterp;
        nCnt = nCnt + 1;
        
        nbPT(nCnt) = I(ii+2, jj);
        nCnt = nCnt + 1;
        
        horInterp1 = I(ii+1, jj-2) + 0.5858*( I(ii+1, jj-1) - I(ii+1, jj-2) );
        horInterp2 = I(ii+2, jj-2) + 0.5858*( I(ii+2, jj-1) - I(ii+2, jj-1) );
        verInterp = horInterp1 + 0.4142*( horInterp2 - horInterp1 );
        nbPT(nCnt) = verInterp;
                
        
        for iCnt = 1:nCnt
            if( nbPT(iCnt) >= I(ii, jj) )
                MatLBP(ii-rad, jj-rad) = MatLBP(ii-rad, jj-rad) + 2^(nCnt-iCnt);
            end
        end
    end
end



% 计算LBP直方图
histLBP = zeros(1, 59); % 对于(8,2)的uniform直方图共有59个收集箱

for ii = 1:m-2*rad
    for jj = 1:n-2*rad
        histLBP( vecLBPMap( MatLBP(ii, jj)+1 ) ) = histLBP( vecLBPMap( MatLBP(ii, jj)+1 ) ) + 1;
    end
end

getLBPHist.m

function histLBP = getLBPHist(I, r, c, nMB)
% 取得I的分区 LBP 直方图
%
% 输入：r,c --- 了分区的数目，r*c个分区
%       nMB --- MB-LBP 中块的大小
%
% 返回值：histLBP ---  连接 I 的各个分块的LBP直方图而形成的代表 I 的LBP复合特征向量

[m n] = size(I);

% 计算分区的大小
mPartitionSize = floor(m / r);
nPartitionSize = floor(n / c);

for ii = 1:r-1
    for jj = 1:c-1
        Sub = I( (ii-1)*mPartitionSize+1:ii*mPartitionSize, (jj-1)*nPartitionSize+1:jj*nPartitionSize );
%        hist{ii}{jj} = getMBLBPFea( Sub, nMB );
        hist{ii}{jj} = getMBLBPFea_33( Sub, nMB );
    end
end


% 处理最后一行和最后一列
clear Sub
for ii = 1:r-1
    Sub = I( (ii-1)*mPartitionSize+1:ii*mPartitionSize, (c-1)*nPartitionSize+1:n );
%    hist{ii}{c} = getMBLBPFea(Sub, nMB);
    hist{ii}{c} = getMBLBPFea_33( Sub, nMB );
end
clear Sub

for jj = 1:c-1
    Sub = I( (r-1)*mPartitionSize+1:m, (jj-1)*nPartitionSize+1:jj*nPartitionSize );
%    hist{r}{jj} = getMBLBPFea(Sub, nMB);
    hist{r}{jj} = getMBLBPFea_33( Sub, nMB );
end
clear Sub

Sub = I((r-1)*mPartitionSize+1:m, (c-1)*nPartitionSize+1:n);
%hist{r}{c} = getMBLBPFea(Sub, nMB);
hist{r}{c} = getMBLBPFea_33( Sub, nMB );


% 连接各个分块的LBP直方图形成复合特征向量
histLBP = zeros(1, 0);
for ii = 1:r
    for jj = 1:c
        histLBP = [histLBP hist{ii}{jj}];
    end
end

makeLBPMap.m

function vecLBPMap = makeLBPMap
% 生成(8,2)临域uniform LBP直方图的映射关系，即将256个灰度值映射到59个收集箱中，
% 所有的非 uniform 放入一个收集箱中

vecLBPMap = zeros(1, 256); %初始化映射表

bits = zeros(1, 8); %8位二进模式串

nCurBin = 1;

for ii = 0:255
    num = ii;
    
    nCnt = 0;
    
    % 获得灰度num的二进制表示bits
    while (num)
        bits(8-nCnt) = mod(num, 2);
        num = floor( num / 2 );
        nCnt = nCnt + 1;
    end
    
    if IsUniform(bits) % 判断bits是不是uniform模式
        vecLBPMap(ii+1) = nCurBin;% 每个uniform模式分配一个收集箱
        nCurBin = nCurBin + 1;
    else
        vecLBPMap(ii+1) = 59;%所有非uniform模式都放入第59号收集箱        
    end
    
end

% 保存映射表
save('Mat/LBPMap.mat', 'vecLBPMap');





function bUni = IsUniform(bits)
% 判断某一个位串模式 bits 是否是 uniform 模式
%
% 输入：bits --- 二进制LBP模式串
%
% 返回值：bUni --- =1，if bits 是uniform模式串；=2，if bits 不是uniform模式串

n = length(bits);

nJmp = 0; % 位跳变数（0->1 or 1->0）
for ii = 1 : (n-1)
    if( bits(ii) ~= bits(ii+1) )
        nJmp = nJmp+1;
    end
end
if bits(n) ~= bits(1)
    nJmp = nJmp+1;
end

if nJmp > 2
    bUni = false;
else
    bUni = true;
end

dist_LBPFea.m

function distVec = dist_LBPFea(LBPFea1, LBPFea2)
% compute the difference vector between 2 LBP features

nLBPFea = 59;
nLen = length(LBPFea1);

nBlock = nLen / nLBPFea;

distVec = zeros(1, nBlock);

for iBlock = 1:nBlock
    % 计算对应块的分块 LBP 直方图距离
    distVec(iBlock) = dist_LBPHist( LBPFea1( (iBlock-1)*nLBPFea+1:iBlock*nLBPFea ), LBPFea2( (iBlock-1)*nLBPFea+1:iBlock*nLBPFea ) );
end

dist_LBPHist.m

function dist = dist_LBPHist(hist1, hist2)
% compute the Chi square statistic distance between 2 LBP histograms

nLBPFea = 59; % nBin


dist = 0;

for iBin = 1:nLBPFea
    dist = dist + (hist1(iBin) - hist2(iBin))^2 / (hist1(iBin)+hist2(iBin)+eps);
%    id dist
end