LBP特征提取原理及代码实现

1 背景

LBP（Local Binary Pattern，局部二值模式）是一种用来描述图像局部纹理特征的算子，具有旋转不变形和灰度值不变形等显著优点。主要用于纹理特征提取，在人脸识别部分有较好的效果。

2 LBP特征原理

2.1概述

　　　　从94年T. Ojala, M.Pietikäinen, 和D. Harwood提出至今，LBP大致经历了三个版本。下面按照时间顺序进行介绍。(下面的First-LBP、Second-LBP..均为自己定义，只是为写技术文档方便，LBP作者并未定义。)

2.2First-LBP原理

　　　　最初的LBP算子通过定义一个3x3的窗口，以窗口内中心点的像素值为标准，对比窗口内另8个点像素值的大小，大于为1，小于为0。8个点形成一个二进制数字（通常转换为十进制表示）即为中心点的LBP特征值。详细计算如下图：

通过上面得到的LBP算子具有很多缺点，后研究人员在LBP基础上进行不断改进。

2.3 Second-LBP原理

原始LBP算子计算区域为像素点的周围8个点，在图像尺寸发生改变时会出现很大的偏差，不能正确反映像素点周围的纹理信息。为适应不同尺寸纹理特征，LBP原作者将圆形邻域代替正方形邻域。同时增加了旋转不变的特性，在对LBP特征值的存储部分，也进行了改进。详细如下文。

2.3.1圆形LBP算子

圆形LBP算子计算中，以像素点为圆心，R为半径，提取半径上P个采样点，根据2.2中像素值比较方法，进行像素值大小的比较，得到该点的LBP特征值。其中提取采样点的方法如下：

x_{t =}x_{d + Rcos(2πp/P)}

　　　　　y_{t =} y_{d - Rsin(2πp/P)}

　　（x_t,y_t）为某个采样点，(x_d,y_d)为邻域中心点，p为第P个采样点，P为采样点的个数。得到采样点的坐标可能为小数，改进后的LBP采用双线性插值法进行计算该点的像素值：

几种不同半径不同采样点数量的LBP算子：

2.3.2旋转不变LBP特征：

　上面通过采取圆形邻域的计算，一定程度上削弱了尺度改变的影响。研究人员在上面的基础上进一步扩展，使具备旋转不变的特征。

　　首先，在确定半径大小和采样点数目后，不断旋转圆形邻域内采样点的位置，得到一系列的LBP特征值，从这些LBP特征值中选择最小的值作为LBP中心像素点的LBP特征值，具体如下图：

通过不断旋转，取最小值，使具备旋转不变特性。

2.3.3 Uniform Pattern LBP特征：

Uniform Pattern LBP特征也称为等价模式或均匀模式。对LBP特征值的存储方式上，进行了优化。详细如下。

假设对于半径为R的圆形邻域内提取P个采样点，会产生2^p种二进制表达方法，随着邻域内采样点数目的增加，二进制模式的种类以指数形式增加，不利于LBP特征值的存储、提取、分类和识别。LBP原作者提出一种“等价模式”对LBP算子进行降维。详细如下。

在实际图像中，绝大多数LBP模式只包括从0到1或从0到1的转变，LBP原作者将“等价模式”定义为当某个LBP特征值所对应的二进制数从0到1或从1到0的转变最多有两次时，该LBP所对应的二进制就称为一个等价模式。如00000000（0次跳变）、00000011（1次跳变）、10001111

（2次跳变）均为等价模式类。除等价模式类外均归为混合模式类。上述算法，使得模式数量由原来的2^p种减少为P（P-1）+2+1种（P代表采样点的数量）。

实例介绍：

如采样点数为8，即256种LBP特征值，根据等价模式可分为59类：跳变0次——2个，跳变1次——0个，跳变2次——56个，跳变3….8次——1个。（跳变1次为0个是因为LBP作者把LBP二进制数字看做一个圆性的序列，故跳变1次为0个）

2.3.4 MB-LBP特征：

MB-LBP特征，全称为Multiscale Block LBP,由中科院的研究人员研究发表，原理与HOG特征提取有相似之处，接扫MB-LBP仅用于了解，下面是原理介绍。

首先将图像分为分为多个块，再将每个小块分成多个区域，每个区域的灰度值为该区域内灰度值的平均值。在一个块内，将中心区域的灰度值大小与周围区域的灰度值大小进行比较形成LBP特征值。如下图：

作者对得到的MB-LBP特征值同样进行均值编码。首先，对得到的特征值采用直方图进行表示，计算每一种特征值的数量，进行排序，将排序在前63为的特征值看作是等价模式类，其他的为混合模式类，共64类。

2.4 Third-LBP原理

LBP的最后一步改进为LBPH即LBP特征统计直方图的使用，可用于机器学习特征的提取。这种表示方法由Ahonen等人提出，将LBP特征图像分成m个局部块，提取每个局部块的直方图，并依次连接在一起形成LBP特征的统计直方图。具体过程如下：

（1）计算图像中每一像素点的LBP特征值。

（2）图像进行分成多块。（Opencv中默认将LBP特征

图像分为8行8列64块区域。

（3）计算每块区域的LBP特征值的直方图，并将直方图进行归一化。（横坐标为LBP特征值的表示方式，纵坐标为数量）

（4）将上面计算的每块区域特征图像的直方图按顺序依次排列成一行，形成LBP特征向量。

（5）用机器学习方法对LBP特征向量进行训练。

举例说明LBPH的维度：
　　　　采样点为8个，如果用的是原始的LBP或Extended LBP特征，其LBP特征值的模式为256种，则一幅图像的LBP特征向量维度为：64*256=16384维，而如果使用的UniformPatternLBP特征，其LBP值的模式为59种，其特征向量维度为：64*59=3776维，可以看出，使用等价模式特征，其特征向量的维度大大减少，这意味着使用机器学习方法进行学习的时间将大大减少，而性能上没有受到很大影响。

代码环节：

　　目前，Opencv内，暂无LBP简便的接口函数，如若使用，可观看下面源代码或使用下面代码。

　　LBP在人脸识别部分的应用：

#include<iostream>
#include<opencv2/opencv.hpp>
#include<opencv2/highgui/highgui.hpp>
#include<opencv2/imgproc/imgproc.hpp>
#include"opencv2/face.hpp"

using namespace std;
using namespace cv;

int main(int argc,char* argv[])
{
    vector<Mat> images;
    vector<int> labels;
    char tdt[10];
    for(int i=1;i<8;i++)
    {
        sprintf(tdt,"0%d.jpg",i);
        Mat image = imread(tdt);
        Mat grayImage;
        cvtColor(image,grayImage,COLOR_BGR2GRAY);
        images.push_back(grayImage);
        labels.push_back(1);
    }
    for(int i=8;i<12;i++)
    {
        sprintf(tdt,"0%d.jpg",i);
        Mat image = imread(tdt);
        Mat grayImage;
        cvtColor(image,grayImage,COLOR_BGR2GRAY);
        images.push_back(grayImage);
        labels.push_back(2);
    }
    Ptr<cv::face::FaceRecognizer> p ;
    p->train(images,labels);
    Mat test= imread("12.jpg");
    Mat grayImage;
    cvtColor(test,grayImage,COLOR_BGR2GRAY);
    int result = p->predict(grayImage);
    cout<<result<<endl;
    system("pause");
    return 0;
}

　　原始LBP特征计算代码：

template <typename _tp>
void getOriginLBPFeature(InputArray _src,OutputArray _dst)
{
    Mat src = _src.getMat();
    _dst.create(src.rows-2,src.cols-2,CV_8UC1);
    Mat dst = _dst.getMat();
    dst.setTo(0);
    for(int i=1;i<src.rows-1;i++)
    {
        for(int j=1;j<src.cols-1;j++)
        {
            _tp center = src.at<_tp>(i,j);
            unsigned char lbpCode = 0;
            lbpCode |= (src.at<_tp>(i-1,j-1) > center) << 7;
            lbpCode |= (src.at<_tp>(i-1,j  ) > center) << 6;
            lbpCode |= (src.at<_tp>(i-1,j+1) > center) << 5;
            lbpCode |= (src.at<_tp>(i  ,j+1) > center) << 4;
            lbpCode |= (src.at<_tp>(i+1,j+1) > center) << 3;
            lbpCode |= (src.at<_tp>(i+1,j  ) > center) << 2;
            lbpCode |= (src.at<_tp>(i+1,j-1) > center) << 1;
            lbpCode |= (src.at<_tp>(i  ,j-1) > center) << 0;
            dst.at<uchar>(i-1,j-1) = lbpCode;
        }
    }
}

　　　　圆形LBP特征计算：

template <typename _tp>
void getCircularLBPFeature(InputArray _src,OutputArray _dst,int radius,int neighbors)
{
    Mat src = _src.getMat();
    //LBP特征图像的行数和列数的计算要准确
    _dst.create(src.rows-2*radius,src.cols-2*radius,CV_8UC1);
    Mat dst = _dst.getMat();
    dst.setTo(0);
    //循环处理每个像素
    for(int i=radius;i<src.rows-radius;i++)
    {
        for(int j=radius;j<src.cols-radius;j++)
        {
            //获得中心像素点的灰度值
            _tp center = src.at<_tp>(i,j);
            unsigned char lbpCode = 0;
            for(int k=0;k<neighbors;k++)
            {
                //根据公式计算第k个采样点的坐标，这个地方可以优化，不必每次都进行计算radius*cos，radius*sin
                float x = i + static_cast<float>(radius *                     cos(2.0 * CV_PI * k / neighbors));
                float y = j - static_cast<float>(radius *                     sin(2.0 * CV_PI * k / neighbors));
                //根据取整结果进行双线性插值，得到第k个采样点的灰度值

                //1.分别对x，y进行上下取整
                int x1 = static_cast<int>(floor(x));
                int x2 = static_cast<int>(ceil(x));
                int y1 = static_cast<int>(floor(y));
                int y2 = static_cast<int>(ceil(y));

                //2.计算四个点(x1,y1),(x1,y2),(x2,y1),(x2,y2)的权重
                //下面的权重计算方式有个问题，如果四个点都相等，则权重全为0，计算出来的插值为0
                //float w1 = (x2-x)*(y2-y); //(x1,y1)
                //float w2 = (x2-x)*(y-y1); //(x1,y2)
                //float w3 = (x-x1)*(y2-y); //(x2,y1)
                //float w4 = (x-x1)*(y-y1); //(x2,y2)

                //将坐标映射到0-1之间
                float tx = x - x1;
                float ty = y - y1;
                //根据0-1之间的x，y的权重计算公式计算权重
                float w1 = (1-tx) * (1-ty);
                float w2 =    tx  * (1-ty);
                float w3 = (1-tx) *    ty;
                float w4 =    tx  *    ty;
                //3.根据双线性插值公式计算第k个采样点的灰度值
                float neighbor = src.at<_tp>(x1,y1) * w1 + src.at<_tp>(x1,y2) *w2                     + src.at<_tp>(x2,y1) * w3 +src.at<_tp>(x2,y2) *w4;
                //通过比较获得LBP值，并按顺序排列起来
                lbpCode |= (neighbor>center) <<(neighbors-k-1);
            }
            dst.at<uchar>(i-radius,j-radius) = lbpCode;
        }
    }
}
//圆形LBP特征计算，效率优化版本，声明时默认neighbors=8
template <typename _tp>
void getCircularLBPFeatureOptimization(InputArray _src,OutputArray _dst,int radius,int neighbors)
{
    Mat src = _src.getMat();
    //LBP特征图像的行数和列数的计算要准确
    _dst.create(src.rows-2*radius,src.cols-2*radius,CV_8UC1);
    Mat dst = _dst.getMat();
    dst.setTo(0);
    for(int k=0;k<neighbors;k++)
    {
        //计算采样点对于中心点坐标的偏移量rx，ry
        float rx = static_cast<float>(radius * cos(2.0 * CV_PI * k / neighbors));
        float ry = -static_cast<float>(radius * sin(2.0 * CV_PI * k / neighbors));
        //为双线性插值做准备
        //对采样点偏移量分别进行上下取整
        int x1 = static_cast<int>(floor(rx));
        int x2 = static_cast<int>(ceil(rx));
        int y1 = static_cast<int>(floor(ry));
        int y2 = static_cast<int>(ceil(ry));
        //将坐标偏移量映射到0-1之间
        float tx = rx - x1;
        float ty = ry - y1;
        //根据0-1之间的x，y的权重计算公式计算权重，权重与坐标具体位置无关，与坐标间的差值有关
        float w1 = (1-tx) * (1-ty);
        float w2 =    tx  * (1-ty);
        float w3 = (1-tx) *    ty;
        float w4 =    tx  *    ty;
        //循环处理每个像素
        for(int i=radius;i<src.rows-radius;i++)
        {
            for(int j=radius;j<src.cols-radius;j++)
            {
                //获得中心像素点的灰度值
                _tp center = src.at<_tp>(i,j);
                //根据双线性插值公式计算第k个采样点的灰度值
                float neighbor = src.at<_tp>(i+x1,j+y1) * w1 + src.at<_tp>(i+x1,j+y2) *w2                     + src.at<_tp>(i+x2,j+y1) * w3 +src.at<_tp>(i+x2,j+y2) *w4;
                //LBP特征图像的每个邻居的LBP值累加，累加通过与操作完成，对应的LBP值通过移位取得
                dst.at<uchar>(i-radius,j-radius) |= (neighbor>center) <<(neighbors-k-1);
            }
        }
    }
}

　　　　旋转不变圆形LBP特征计算：

template <typename _tp>
void getRotationInvariantLBPFeature(InputArray _src,OutputArray _dst,int radius,int neighbors)
{
    Mat src = _src.getMat();
    //LBP特征图像的行数和列数的计算要准确
    _dst.create(src.rows-2*radius,src.cols-2*radius,CV_8UC1);
    Mat dst = _dst.getMat();
    dst.setTo(0);
    for(int k=0;k<neighbors;k++)
    {
        //计算采样点对于中心点坐标的偏移量rx，ry
        float rx = static_cast<float>(radius * cos(2.0 * CV_PI * k / neighbors));
        float ry = -static_cast<float>(radius * sin(2.0 * CV_PI * k / neighbors));
        //为双线性插值做准备
        //对采样点偏移量分别进行上下取整
        int x1 = static_cast<int>(floor(rx));
        int x2 = static_cast<int>(ceil(rx));
        int y1 = static_cast<int>(floor(ry));
        int y2 = static_cast<int>(ceil(ry));
        //将坐标偏移量映射到0-1之间
        float tx = rx - x1;
        float ty = ry - y1;
        //根据0-1之间的x，y的权重计算公式计算权重，权重与坐标具体位置无关，与坐标间的差值有关
        float w1 = (1-tx) * (1-ty);
        float w2 =    tx  * (1-ty);
        float w3 = (1-tx) *    ty;
        float w4 =    tx  *    ty;
        //循环处理每个像素
        for(int i=radius;i<src.rows-radius;i++)
        {
            for(int j=radius;j<src.cols-radius;j++)
            {
                //获得中心像素点的灰度值
                _tp center = src.at<_tp>(i,j);
                //根据双线性插值公式计算第k个采样点的灰度值
                float neighbor = src.at<_tp>(i+x1,j+y1) * w1 + src.at<_tp>(i+x1,j+y2) *w2                     + src.at<_tp>(i+x2,j+y1) * w3 +src.at<_tp>(i+x2,j+y2) *w4;
                //LBP特征图像的每个邻居的LBP值累加，累加通过与操作完成，对应的LBP值通过移位取得
                dst.at<uchar>(i-radius,j-radius) |= (neighbor>center) <<(neighbors-k-1);
            }
        }
    }
    //进行旋转不变处理
    for(int i=0;i<dst.rows;i++)
    {
        for(int j=0;j<dst.cols;j++)
        {
            unsigned char currentValue = dst.at<uchar>(i,j);
            unsigned char minValue = currentValue;
            for(int k=1;k<neighbors;k++)
            {
    //循环左移
                unsigned char temp = (currentValue>>(neighbors-k)) | (currentValue<<k);
                if(temp < minValue)
                {
                    minValue = temp;
                }
            }
            dst.at<uchar>(i,j) = minValue;
        }
    }
}

　　　　均匀（等价）模式：

template <typename _tp>
void getUniformPatternLBPFeature(InputArray _src,OutputArray _dst,int radius,int neighbors)
{
    Mat src = _src.getMat();
    //LBP特征图像的行数和列数的计算要准确
    _dst.create(src.rows-2*radius,src.cols-2*radius,CV_8UC1);
    Mat dst = _dst.getMat();
    dst.setTo(0);
    //LBP特征值对应图像灰度编码表，直接默认采样点为8位
    uchar temp = 1;
    uchar table[256] = {0};
    for(int i=0;i<256;i++)
    {
        if(getHopTimes(i)<3)
        {
            table[i] = temp;
            temp++;
        }
    }
    //是否进行UniformPattern编码的标志
    bool flag = false;
    //计算LBP特征图
    for(int k=0;k<neighbors;k++)
    {
        if(k==neighbors-1)
        {
            flag = true;
        }
        //计算采样点对于中心点坐标的偏移量rx，ry
        float rx = static_cast<float>(radius * cos(2.0 * CV_PI * k / neighbors));
        float ry = -static_cast<float>(radius * sin(2.0 * CV_PI * k / neighbors));
        //为双线性插值做准备
        //对采样点偏移量分别进行上下取整
        int x1 = static_cast<int>(floor(rx));
        int x2 = static_cast<int>(ceil(rx));
        int y1 = static_cast<int>(floor(ry));
        int y2 = static_cast<int>(ceil(ry));
        //将坐标偏移量映射到0-1之间
        float tx = rx - x1;
        float ty = ry - y1;
        //根据0-1之间的x，y的权重计算公式计算权重，权重与坐标具体位置无关，与坐标间的差值有关
        float w1 = (1-tx) * (1-ty);
        float w2 =    tx  * (1-ty);
        float w3 = (1-tx) *    ty;
        float w4 =    tx  *    ty;
        //循环处理每个像素
        for(int i=radius;i<src.rows-radius;i++)
        {
            for(int j=radius;j<src.cols-radius;j++)
            {
                //获得中心像素点的灰度值
                _tp center = src.at<_tp>(i,j);
                //根据双线性插值公式计算第k个采样点的灰度值
                float neighbor = src.at<_tp>(i+x1,j+y1) * w1 + src.at<_tp>(i+x1,j+y2) *w2                     + src.at<_tp>(i+x2,j+y1) * w3 +src.at<_tp>(i+x2,j+y2) *w4;
                //LBP特征图像的每个邻居的LBP值累加，累加通过与操作完成，对应的LBP值通过移位取得
                dst.at<uchar>(i-radius,j-radius) |= (neighbor>center) <<(neighbors-k-1);
                //进行LBP特征的UniformPattern编码
                if(flag)
                {
                    dst.at<uchar>(i-radius,j-radius) = table[dst.at<uchar>(i-radius,j-radius)];
                }
            }
        }
    }
}
//计算跳变次数
int getHopTimes(int n)
{
    int count = 0;
    bitset<8> binaryCode = n;
    for(int i=0;i<8;i++)
    {
        if(binaryCode[i] != binaryCode[(i+1)%8])
        {
            count++;
        }
    }
    return count;
}

MB-LBP代码：

//MB-LBP特征的计算
void getMultiScaleBlockLBPFeature(InputArray _src,OutputArray _dst,int scale)
{
    Mat src = _src.getMat();
    Mat dst = _dst.getMat();
    //定义并计算积分图像
    int cellSize = scale / 3;
    int offset = cellSize / 2;
    Mat cellImage(src.rows-2*offset,src.cols-2*offset,CV_8UC1);
    for(int i=offset;i<src.rows-offset;i++)
    {
        for(int j=offset;j<src.cols-offset;j++)
        {
            int temp = 0;
            for(int m=-offset;m<offset+1;m++)
            {
                for(int n=-offset;n<offset+1;n++)
                {
                    temp += src.at<uchar>(i+n,j+m);
                }
            }
            temp /= (cellSize*cellSize);
            cellImage.at<uchar>(i-cellSize/2,j-cellSize/2) = uchar(temp); 
        }
    }
    getOriginLBPFeature<uchar>(cellImage,dst);
}

MB-LBP编码的计算：

//求SEMB-LBP
void SEMB_LBPFeature(InputArray _src,OutputArray _dst,int scale)
{
    Mat dst=_dst.getMat();
    Mat MB_LBPImage;
    getMultiScaleBlockLBPFeature(_src,MB_LBPImage,scale);
    //imshow("dst",dst);
    Mat histMat;
    int histSize = 256;
    float range[] = {float(0),float(255)};
    const float* ranges = {range};
    //计算LBP特征值0-255的直方图
    calcHist(&MB_LBPImage,1,0,Mat(),histMat,1,&histSize,&ranges,true,false);
    histMat.reshape(1,1);
    vector<float> histVector(histMat.rows*histMat.cols);
    uchar table[256];
    memset(table,64,256);
    if(histMat.isContinuous())
    {
        //histVector = (int *)(histMat.data);
        //将直方图histMat变为vector向量histVector
        histVector.assign((float*)histMat.datastart,(float*)histMat.dataend);
        vector<float> histVectorCopy(histVector);
        //对histVector进行排序，即对LBP特征值的数量进行排序，降序排列
        sort(histVector.begin(),histVector.end(),greater<float>());
        for(int i=0;i<63;i++)
        {
            for(int j=0;j<histVectorCopy.size();j++)
            {
                if(histVectorCopy[j]==histVector[i])
                {
                    //得到类似于Uniform的编码表
                    table[j]=i;
                }
            }
        }
    }
    dst = MB_LBPImage;
    //根据编码表得到SEMB-LBP
    for(int i=0;i<dst.rows;i++)
    {
        for(int j=0;j<dst.cols;j++)
        {
            dst.at<uchar>(i,j) = table[dst.at<uchar>(i,j)];
        }
    }
}

LBPH特征图像的计算：

//计算LBP特征图像的直方图LBPH
Mat getLBPH(InputArray _src,int numPatterns,int grid_x,int grid_y,bool normed)
{
    Mat src = _src.getMat();
    int width = src.cols / grid_x;
    int height = src.rows / grid_y;
    //定义LBPH的行和列，grid_x*grid_y表示将图像分割成这么些块，numPatterns表示LBP值的模式种类
    Mat result = Mat::zeros(grid_x * grid_y,numPatterns,CV_32FC1);
    if(src.empty())
    {
        return result.reshape(1,1);
    }
    int resultRowIndex = 0;
    //对图像进行分割，分割成grid_x*grid_y块，grid_x，grid_y默认为8
    for(int i=0;i<grid_x;i++)
    {
        for(int j=0;j<grid_y;j++)
        {
            //图像分块
            Mat src_cell = Mat(src,Range(i*height,(i+1)*height),Range(j*width,(j+1)*width));
            //计算直方图
            Mat hist_cell = getLocalRegionLBPH(src_cell,0,(numPattern-1),true);
            //将直方图放到result中
            Mat rowResult = result.row(resultRowIndex);
            hist_cell.reshape(1,1).convertTo(rowResult,CV_32FC1);
            resultRowIndex++;
        }
    }
    return result.reshape(1,1);
}
//计算一个LBP特征图像块的直方图
Mat getLocalRegionLBPH(const Mat& src,int minValue,int maxValue,bool normed)
{
    //定义存储直方图的矩阵
    Mat result;
    //计算得到直方图bin的数目，直方图数组的大小
    int histSize = maxValue - minValue + 1;
    //定义直方图每一维的bin的变化范围
    float range[] = { static_cast<float>(minValue),static_cast<float>(maxValue + 1) };
    //定义直方图所有bin的变化范围
    const float* ranges = { range };
    //计算直方图，src是要计算直方图的图像，1是要计算直方图的图像数目，0是计算直方图所用的图像的通道序号，从0索引
    //Mat()是要用的掩模，result为输出的直方图，1为输出的直方图的维度，histSize直方图在每一维的变化范围
    //ranges，所有直方图的变化范围（起点和终点）
    calcHist(&src,1,0,Mat(),result,1,&histSize,&ranges,true,false);
    //归一化
    if(normed)
    {
        result /= (int)src.total();
    }
    //结果表示成只有1行的矩阵
    return result.reshape(1,1);
}

posted @ 2020-03-06 14:14 无趣的鱼阅读(9380) 评论(0) 编辑收藏举报

刷新页面返回顶部

BLUEFISH