最近在一些机器视觉群中的一些小伙伴们多次问到opencv是否集成了LBP算法，据我了解opencv没有单独的LBP特征描述算法实现，都是和一些应用结合，如人脸识别，检测等，这些都是一些论文的研究成果，针对于特定的应用，这对于想将LBP特征描述用到自己的应用中的伙伴来说，或许不太方便。opencv也没有一个单独的特征描述这样的一个module，这或许限制了使用opencv的灵活性，而且个人体会视觉问题最终根本的问题都落在了特征描述这样的一个最基本的问题，包括近几年很火的深度学习，其解决的最根本问题也是特征表达。我们已有的机器学习理论非常的丰富，一个区分能力强，不变性强的特征往往比改进机器学习算法的到的识别效果要好的多。期待opencv开源者在特征描述这块增加更多的内容，也期待自己某天也作出一点贡献。不过现阶段的特征描述与我们语义上的特征在实际应用中差别还很大，多数都是更具经验和知识进行手工描述的，近些年也有一些data driven的方法，但是还有很长的一段路要走,还需要大量研究者作出新的贡献，或许某天我们能够将人的视觉系统如何处理和存储信息一步一步的描述清楚，机器视觉问题才能最终得到解决。呵呵，这都是自己的一些浅显看法，希望大家批评指正。下面开始进入正题吧。
局部二值模式（local binary pattern,LBP）算子是一个描述图像局部空间结构的非参数模型算子,是由T.Ojala于1996年首先提出的,而且Ojala等又于2002年对LBP算子在纹理分类中的高区分能力予以了证明。LBP算法思想简单容易理解、计算复杂度小、对单调的灰度变换不敏感并且能够很好地描述图像的局部纹理特征，广泛的应用到识别领域。LBP算子经过了一系列的改进，到目前为止，或许CLBP中提出的3种差异描述算子组合的改进是最好的。这篇blog在CLBP基础之上又提出了4种其它的差异算子，其中有些算子的组合也作为研究成果发表了，这里基于OpenCV2将7种算子用C++实现，共享出来，你可以用你的“想象力”将这些算子组合作为特征提取，应用到实际应用中。下面先对LBP及相关改进进行简单描述，然后给出7种差异算子的源码及实验结果。
1.基本LBP
下图给出了一个基本LBP算子，应用LBP算子的过程类似于滤波过程中的模板操作。逐行扫描图像中的每一个像素点，以改点的灰度值作为阈值，对其周围3×3的8邻域进行二值化，按照一定的顺序将二值化的结果组成一个8位二进制数，以此二进制数的值作为改点的响应。
例如对与上图中的3×3区域的中心点，以其灰度值88作为阈值，对其8邻域进行二值化，并且从左上角点按照顺时针方向将二值化的结果组成一个二进制数10001011，即十进制的139，作为中心点的响应。在整个逐行扫描过程结束后，会得到一个LBP响应图像，这个响应图像的直方图被称为LBP统计直方图，或LBP直方图，它常常作为后续识别工作的特征，以此也叫做LBP特征。
LBP的主要思想是以某一点与其邻域像素的相对灰度作为响应，正是这种相对机制使LBP算子对于单调的灰度变化具有不变性（对光照变化具有较强的鲁棒性，可用实验结果验证）。
2.圆形邻域的LBP
采用圆形的邻域并结合双线性插值运算使我们可以获得任意半径和任意数目的邻域像素点。例如一个半径为2的8邻域像素的圆形邻域，每一个方格对应一个像素，对于正好处于方格中心的邻域点（左、上、右、下四个黑点），直接以改点所在方格的像素值作为它的值；对于不在像素中心位置的邻域点，通过双线性插值确定其值。如下图所示：
3.均匀LBP
由于LBP直方图大多都是针对图像中的各个分区分别计算的，对于一个普通大小的分块区域，标准LBP算子得到的二值模式数目较多，而实际位于该分块区域中的像素数目却相对较少，这将会得到一个过于稀疏的直方图，从而使直方图失去统计意义。以此应设法减少一些冗余的LBP模式，同时又保留足够的具有重要描绘能力的模式。
正是基于上述考虑，研究者提出了均匀模式（uniform patterns）的概念，这是对LBP算子的由一重大改进。对于一个局部二值模式，在将其二进制位串视为循环的情况下，如果其中包含的从0到1或者从1到0转变不多于两个，则称这个局部二值模式为均匀模式。例如，模式00000000（0个转变），01110000（2个转变）和11001111（2个转变）都是均匀模式，而模式11001001（4个转变）和01010011（6个转变）则不是。
均匀模式的意义在于：在随后的LBP直方图的计算过程中，只为均匀模式分配独立的直方图收集箱，而所有的非均匀模式都放入一个公用收集箱，这就使得LBP特征的数目大大减少。一般来说，保留的均匀模式往往都是反映重要信息的那些模式，而那些非均匀模式中过多的转变往往有噪声引起，不具有良好的统计意义。
假设图像分块区域大小为18×20，像素总数为360.如果采用8邻域像素的标准LBP算子，收集箱数目为256个，平均到每个收集箱的像素数目还不到2个，没有统计意义；而均匀LBP算子的收集箱数目为59（58个均匀模式收集箱加上1个非均匀模式收集箱），平均到每个收集箱中的像素个数为6个左右，更具有统计意义。对16邻域像素而言，标准LBP算子和均匀LBP算子的收集箱数目分别为65536和243。
4.旋转不变的LBP
不断旋转圆形邻域得到一系列初始定义的LBP值，取其最小值作为该邻域的LBP值。通过引入旋转不变的定义，使得LBP算子对于图像旋转表现得更为鲁棒，并且LBP统计模式的数量进一步减少。
5.旋转不变的均匀LBP
将旋转不变的均匀LBP，进一步减少了LBP统计模式的数量。
6.完整LBP（CompleteLBP,CLBP）
CLBP是12年被提出的，用三种局部差异算子描述纹理特征：局部灰度差异描述算子（CLBP-Sign, CLBP_S）、局部梯度差异描述算子（CLBP-Magnitude, CLBP_M）以及中心像素点描述算子（CLBP-Center, CLBP_C）。三个算子的计算方式如下：
其中，gc为中心像素灰度值，gp为中心像素点邻域灰度值，N为图像局部中心像素个数。从公式中可以看出：CLBP_SP,R即为传统意义上的LBP；CLBP_MP,R通过两像素点的灰度差异幅值与全局灰度差异幅值的均值比较，描述了局部梯度差异信息；CLBP_CP,R反应中心像素点的灰度信息。相比于传统LBP及其变种，此三种描述算子联合构成的CLBP_SMC对纹理的描述更加精细，对纹理的识别准确率有了大幅度提高。
7.基于OpenCV2实现7种差异描述算子
在CLBP中的3种算子基础上提出了另外4种差异描述算子，它们分别在CLBP_M和CLBP_C算子基础之上各提出了另外2种类CLBP_M和CLBP_C差异描述算子。如果你在熟悉CLBP的基础之上理解这些算子是相对较容易的，可以参考相关论文：《A completed modeling of local binary pattern operator for textureclassification》，《用于纹理特征提取的改进的LBP算法》。下面就直接给出代码好了。
class ICLBP {
private:
double global_average_difference(const cv::Mat& src, int radius=1, int neighbors=8) {
assert(radius>0 && neighbors>0);
//compute global average difference
double averageDiff = 0;
for(int n=0; n<neighbors; n++){
// position of sample point
float x = static_cast<float>(-radius * sin(2.0*CV_PI*n/static_cast<float>(neighbors)));
float y = static_cast<float>(radius * cos(2.0*CV_PI*n/static_cast<float>(neighbors)));
// flooring and ceiling
int fx = static_cast<int>(floor(x));
int fy = static_cast<int>(floor(y));
int cx = static_cast<int>(ceil(x));
int cy = static_cast<int>(ceil(y));
// decimals
float ty = y - fy;
float tx = x - fx;
// weight of interpolation
float w1 = (1 - tx) * (1 - ty);
float w2 = tx * (1 - ty);
float w3 = (1 - tx) * ty;
float w4 = tx * ty;
// processing image data
for(int i=radius; i < src.rows-radius;i++){
for(int j=radius;j < src.cols-radius;j++) {
// compute interpolation value
float t = static_cast<float>(w1*src.at<uchar>(i+fy,j+fx) + w2*src.at<uchar>(i+fy,j+cx)
+ w3*src.at<uchar>(i+cy,j+fx) + w4*src.at<uchar>(i+cy,j+cx));
averageDiff += abs(t-src.at<uchar>(i,j));
}
}//processing data for
}//outer for
averageDiff /= (neighbors*(src.rows-radius)*(src.cols-radius));
return averageDiff;
}
double local_average_difference(const cv::Mat& src, int i, int j, int radius=1, int neighbors=8) {
assert(i>=0 && i<src.cols && j>=0 && j<src.rows);
assert(radius>0 && neighbors>0);
double averageDiff = 0;
for(int n=0; n<neighbors; n++){
// compute position of sample point
float x = static_cast<float>(-radius * sin(2.0*CV_PI*n/static_cast<float>(neighbors)));
float y = static_cast<float>(radius * cos(2.0*CV_PI*n/static_cast<float>(neighbors)));
// flooring and ceiling
int fx = static_cast<int>(floor(x));
int fy = static_cast<int>(floor(y));
int cx = static_cast<int>(ceil(x));
int cy = static_cast<int>(ceil(y));
// decimals
float ty = y - fy;
float tx = x - fx;
// weight of interpolation
float w1 = (1 - tx) * (1 - ty);
float w2 = tx * (1 - ty);
float w3 = (1 - tx) * ty;
float w4 = tx * ty;
// sum of interpolation value
float t = static_cast<float>(w1*src.at<uchar>(i+fy,j+fx) + w2*src.at<uchar>(i+fy,j+cx)
+ w3*src.at<uchar>(i+cy,j+fx) + w4*src.at<uchar>(i+cy,j+cx));
averageDiff += abs(t-src.at<uchar>(i,j));
}
averageDiff /= neighbors;
return averageDiff;
}
double global_average_gray(const cv::Mat& src) {
double whole = 0;
for(int i=0; i<src.rows; i++) {
for(int j=0; j<src.cols; j++) {
whole += src.at<uchar>(i,j);
}
}
whole /= (src.rows)*(src.cols);
return whole;
}
double local_average_gray(const cv::Mat& src, int i, int j, int radius=1, int neighbors=8) {
assert(i>=0 && i<src.cols && j>=0 && j<src.rows);
assert(radius>0 && neighbors>0);
double t = 0;
for(int n=0; n<neighbors; n++){
// compute position of sample point
float x = static_cast<float>(-radius * sin(2.0*CV_PI*n/static_cast<float>(neighbors)));
float y = static_cast<float>(radius * cos(2.0*CV_PI*n/static_cast<float>(neighbors)));
// flooring ceiling
int fx = static_cast<int>(floor(x));
int fy = static_cast<int>(floor(y));
int cx = static_cast<int>(ceil(x));
int cy = static_cast<int>(ceil(y));
// decimals
float ty = y - fy;
float tx = x - fx;
// weight of interpolation
float w1 = (1 - tx) * (1 - ty);
float w2 = tx * (1 - ty);
float w3 = (1 - tx) * ty;
float w4 = tx * ty;
// sum of interpolation value
t += static_cast<float>(w1*src.at<uchar>(i+fy,j+fx) + w2*src.at<uchar>(i+fy,j+cx)
+ w3*src.at<uchar>(i+cy,j+fx) + w4*src.at<uchar>(i+cy,j+cx));
}
t /= neighbors;
return t;
}
public:
//local gray difference code
cv::Mat lbps(const cv::Mat& src, int radius=1, int neighbors=8){
assert(src.data);
assert(radius > 0 && neighbors > 0);
assert((src.type()==CV_8UC1));
assert(radius>0 && neighbors>0);
cv::Mat dst(src.rows, src.cols, CV_8UC1);
dst.setTo(0);
for(int n=0; n<neighbors; n++){
// position of sample point
float x = static_cast<float>(-radius * sin(2.0*CV_PI*n/static_cast<float>(neighbors)));
float y = static_cast<float>(radius * cos(2.0*CV_PI*n/static_cast<float>(neighbors)));
// flooring and ceiling
int fx = static_cast<int>(floor(x));
int fy = static_cast<int>(floor(y));
int cx = static_cast<int>(ceil(x));
int cy = static_cast<int>(ceil(y));
// decimals
float ty = y - fy;
float tx = x - fx;
// weight of interpolation
float w1 = (1 - tx) * (1 - ty);
float w2 = tx * (1 - ty);
float w3 = (1 - tx) * ty;
float w4 = tx * ty;
// processing image data
for(int i=radius; i < src.rows-radius;i++){
for(int j=radius;j < src.cols-radius;j++) {
// compute interpolation value
float t = static_cast<float>(w1*src.at<uchar>(i+fy,j+fx) + w2*src.at<uchar>(i+fy,j+cx)
+ w3*src.at<uchar>(i+cy,j+fx) + w4*src.at<uchar>(i+cy,j+cx));
// coding
dst.at<uchar>(i-radius,j-radius) += ((t > src.at<uchar>(i,j)) ||
(std::abs(t-src.at<uchar>(i,j)) < std::numeric_limits<float>::epsilon())) << n;
}
}//processing data for
}//outer for
return dst;
}
//differnece between local differnce and global average differnce code
cv::Mat lbpm(const cv::Mat& src, int radius=1, int neighbors=8) {
assert(src.data);
assert(radius > 0 && neighbors > 0);
assert((src.type()==CV_8UC1));
assert(radius>0 && neighbors>0);
cv::Mat dst(src.rows, src.cols, CV_8UC1);
dst.setTo(0);
double averageDiff = global_average_difference(src, radius, neighbors);
for(int n=0; n<neighbors; n++){
// position of sample point
float x = static_cast<float>(-radius * sin(2.0*CV_PI*n/static_cast<float>(neighbors)));
float y = static_cast<float>(radius * cos(2.0*CV_PI*n/static_cast<float>(neighbors)));
// flooring and ceiling
int fx = static_cast<int>(floor(x));
int fy = static_cast<int>(floor(y));
int cx = static_cast<int>(ceil(x));
int cy = static_cast<int>(ceil(y));
// decimals
float ty = y - fy;
float tx = x - fx;
// weight of interpolation
float w1 = (1 - tx) * (1 - ty);
float w2 = tx * (1 - ty);
float w3 = (1 - tx) * ty;
float w4 = tx * ty;
// processing image data
for(int i=radius; i < src.rows-radius;i++){
for(int j=radius;j < src.cols-radius;j++) {
// compute interpolation value
float t = static_cast<float>(w1*src.at<uchar>(i+fy,j+fx) + w2*src.at<uchar>(i+fy,j+cx)
+ w3*src.at<uchar>(i+cy,j+fx) + w4*src.at<uchar>(i+cy,j+cx));
// coding
t -= src.at<uchar>(i,j);
dst.at<uchar>(i-radius,j-radius) += ((abs(t) > averageDiff) ||
(std::abs(abs(t)-averageDiff) < std::numeric_limits<float>::epsilon())) << n;
}
}//processing data for
}//outer for
return dst;
}
//differnece between local differnce and local average differnce code
cv::Mat lbpm_local(const cv::Mat& src, int radius=1, int neighbors=8) {
assert(src.data);
assert(radius > 0 && neighbors > 0);
assert((src.type()==CV_8UC1));
assert(radius>0 && neighbors>0);
cv::Mat dst(src.rows, src.cols, CV_8UC1);
dst.setTo(0);
for(int i=radius; i < src.rows-radius;i++){
for(int j=radius;j < src.cols-radius;j++) {
double averageDiff = local_average_difference(src, i, j, radius, neighbors);
for(int n=0; n<neighbors; n++){
// compute position of sample point
float x = static_cast<float>(-radius * sin(2.0*CV_PI*n/static_cast<float>(neighbors)));
float y = static_cast<float>(radius * cos(2.0*CV_PI*n/static_cast<float>(neighbors)));
// flooring and ceiling
int fx = static_cast<int>(floor(x));
int fy = static_cast<int>(floor(y));
int cx = static_cast<int>(ceil(x));
int cy = static_cast<int>(ceil(y));
// decimals
float ty = y - fy;
float tx = x - fx;
// weight of interpolation
float w1 = (1 - tx) * (1 - ty);
float w2 = tx * (1 - ty);
float w3 = (1 - tx) * ty;
float w4 = tx * ty;
// sum of interpolation value
float t = static_cast<float>(w1*src.at<uchar>(i+fy,j+fx) + w2*src.at<uchar>(i+fy,j+cx)
+ w3*src.at<uchar>(i+cy,j+fx) + w4*src.at<uchar>(i+cy,j+cx));
// coding
t -= src.at<uchar>(i,j);
dst.at<uchar>(i-radius,j-radius) += ((abs(t) > averageDiff) ||
(std::abs(abs(t)-averageDiff) < std::numeric_limits<float>::epsilon())) << n;
}
}//for
}//for
return dst;
}
//differnece between local average differnce and local average differnce
cv::Mat lbpm_local_whole(const cv::Mat& src, int radius=1, int neighbors=8) {
assert(src.data);
assert(radius > 0 && neighbors > 0);
assert((src.type()==CV_8UC1));
assert(radius>0 && neighbors>0);
cv::Mat dst(src.rows, src.cols, CV_8UC1);
dst.setTo(0);
double global = global_average_difference(src, radius, neighbors);
for(int i=radius; i < src.rows-radius;i++){
for(int j=radius;j < src.cols-radius;j++) {
double local = local_average_difference(src, i, j, radius, neighbors);
if(local > global) {
dst.at<uchar>(i, j) = 255;
}else {
dst.at<uchar>(i, j) = 0;
}
}
}
return dst;
}
//diffrence between local gray and global average gray
cv::Mat lbpc(const cv::Mat& src) {
assert(src.data);
assert((src.type()==CV_8UC1));
cv::Mat dst(src.rows, src.cols, CV_8UC1);
dst.setTo(0);
double whole = global_average_gray(src);
for(int i=0; i<src.rows; i++) {
for(int j=0; j<src.cols; j++) {
if(src.at<uchar>(i,j) > whole) {
dst.at<uchar>(i, j) = 255;
}else {
dst.at<uchar>(i, j) = 0;
}
}
}
return dst;
}
//difference between local gray and local average gray
cv::Mat lbpc_local(const cv::Mat& src, int radius=1, int neighbors=8) {
assert(src.data);
assert(radius > 0 && neighbors > 0);
assert((src.type()==CV_8UC1));
cv::Mat dst(src.rows, src.cols, CV_8UC1);
dst.setTo(0);
for(int i=radius; i < src.rows-radius;i++){
for(int j=radius;j < src.cols-radius;j++) {
double local = local_average_gray(src, i, j, radius, neighbors);
if(src.at<uchar>(i, j) > local) {
dst.at<uchar>(i-radius, j-radius) = 255;
}else {
dst.at<uchar>(i-radius, j-radius) = 0;
}
}
}
return dst;
}
//difference between local average gray and local average gray
cv::Mat lbpc_local_whole(const cv::Mat& src, int radius=1, int neighbors=8) {
assert(src.data);
assert(radius > 0 && neighbors > 0);
assert((src.type()==CV_8UC1));
assert(radius>0 && neighbors>0);
cv::Mat dst(src.rows, src.cols, CV_8UC1);
dst.setTo(0);
double whole = global_average_gray(src);
for(int i=radius; i < src.rows-radius;i++){
for(int j=radius;j < src.cols-radius;j++) {
double local = local_average_gray(src, i, j, radius, neighbors);
if(local > whole) {
dst.at<uchar>(i-radius, j-radius) = 255;
}else {
dst.at<uchar>(i-radius, j-radius) = 0;
}
}
}
return dst;
}
};
8.LBP特征描述实验结果
lbps:
lbpm:
lbpm_local:
lbpm_local_whole:
lbpc:
lbpc_local:
lbpc_local_whole:
9.LBP统计特征
为了方便统计均匀LBP特征和旋转不变LBP特征及旋转不变均匀LBP特征模式，先将常规的LBP特征模式映射到相应均匀，旋转不变及旋转不变均模式，方便后续的特征组合的统计特征的提取工作。下面的mapping结构和相关函数完成了这个工作。
struct Mapping {
int* table;
int samples;
int num;
Mapping() {
samples =0;
num = 0;
table = NULL;
}
~Mapping() {
if (table != NULL)
delete[] table;
}
private:
Mapping& operator=(const Mapping&);
Mapping(const Mapping&);
};
static int GetTransU2(int i, int samples) {
int trans = 0;
int temp = i;
temp = temp>>1;
int pre = i - (temp<<1);
int next;
int j = 2;
do {
int temp2 = temp;
temp2 = temp2>>1;
next = temp -(temp2<<1);
if(pre != next) {
trans++;
}
pre = next;
temp = temp2;
j++;
}while(j <= samples);
return trans;
}
static int GetMinRi(int i, int samples) {
int rm = i;
int r = i;
for(int j=0; j<samples-1; j++) {
int temp = r;
temp = temp>>1;
int out = r - (temp<<1);
out = out<<(samples-1);
r = out + temp;
if(r < rm) {
rm = r;
}
}
return rm;
}
static int GetNumOf1(int i, int samples) {
int num = 0;
for(int j=0; j<samples; j++) {
int temp = i;
temp = temp>>1;
int out = i - (temp<<1);
num += out;
i = temp;
}
return num;
}
void GetMapping(Mapping& map, int samples=8, std::string mappingType="u2") {
assert(samples >0);
assert(mappingType == "u2" || mappingType == "ri" || mappingType == "riu2");
int newMax = 0; //number of new patterns
int tableNum = 1;
for(int i=0; i<samples; i++) {
tableNum <<= 1;
}
map.table = new int[tableNum];
if(mappingType == "u2") {
newMax = samples*(samples-1) + 3;
int index = 0;
for(int i=0; i<tableNum; i++) {
int trans = GetTransU2(i, samples);
if(trans <= 2) {
map.table[i] = index;
index++;
}else {
map.table[i] = newMax-1;
}
}//for
}else if(mappingType == "ri") {
int *tempTable = new int[tableNum];
for(int i=0; i<tableNum; i++) {
tempTable[i] = -1;
}
for(int i=0; i<tableNum; i++) {
int rm = GetMinRi(i, samples);
if (tempTable[rm] < 0) {
tempTable[rm] = newMax;
newMax = newMax + 1;
}
map.table[i] = tempTable[rm];
}//for
delete[] tempTable;
}else if(mappingType == "riu2") {
newMax = samples + 2;
for(int i=0; i<tableNum; i++) {
int trans = GetTransU2(i, samples);
if(trans <= 2) {
map.table[i] = GetNumOf1(i, samples);
}else {
map.table[i] = samples + 1;
}
}
}
map.samples = samples;
map.num = newMax;
}
10.一个特征组合的统计特征提取例子
//get statistical feature form combination of lbps_r and lbpc_r
static std::vector<double> GetLBP_S_C(const cv::Mat& src, Mapping& Map,
int radius=1, int neighbors=8,
std::string histType="hs") {
assert(src.data != NULL);
assert(radius > 0 && neighbors > 0);
assert((src.type()==CV_8UC1));
assert(histType == "u2" || histType == "ri" || histType == "riu2");
//initialize feature vector
std::vector<double> hist;
for(int i=0; i<Map.num*2; i++) {
hist.push_back(0);
}
ICLBP clbp;
cv::Mat lbps = clbp.lbps(src, radius, neighbors);
cv::Mat lbpc = clbp.lbpc(src, radius, neighbors);
for(int i=1; i<src.rows-1; i++) {
for(int j=1; j<src.cols-1; j++) {
int pos = Map.table[lbps.at<uchar>(i-1, j-1)];
if(lbpc.at<uchar>(i-1, j-1) != 0) {
hist[2*pos] ++;
}else {
hist[2*pos+1]++;
}
}
}
if(histType == "hs") {
double whole = 0;
for(int i=0; i<Map.num*2; i++) {
whole += hist[i];
}
for(int i=0; i<Map.num*2; i++) {
hist[i] /= whole;
}
}
return hist;
}
reference:
《A completed modeling of local binary pattern operator for textureclassification》
《用于纹理特征提取的改进的LBP算法》
时间： 2024-10-14 03:01:52