LBP 源码分析

一、前言：

  因为最近需要实现一个图片分类器，所以又去看了点图像特征提取方面资料。比较常用的图像特征有Harr角点检测，HOG（梯度方向直方图），LBP等等。而我本次采用的是LBP特征。

  LBP优点是相对于其他图像特征来说原理较为简单，而缺点之一是较低版本的matlab和opencv中都没有将其封装为函数。matlab2015b加入了extractLBPFeature函数来提取图像的lbp特征。可是问题1. matlab图像处理的相关代码应该都不能生成c或者c++代码。2.在于提取后的特征向量和按照论文【1】所得到的结果不一致（不知道是我对论文理解有问题还是其他原因）。

  关于LBP特征讲解的博文有以下推荐：

http://blog.csdn.net/stdcoutzyx/article/details/37317863

http://blog.csdn.net/xidianzhimeng/article/details/19634573

二、源码分析：

以下的lbp代码来自：http://www.bytefish.de/blog/local_binary_patterns/， 在他的主页上有对lbp特征的简要介绍。主要的lbp特征实现代码在lbp.cpp和lbpHistogram.cpp两个文件中

2.1:lbp.cpp

#include "lbp.hpp"using namespace cv;//Original LBP template <typename _Tp>void lbp::OLBP_(const Mat& src, Mat& dst) {//设置结果矩阵的行列数dst = Mat::zeros(src.rows - 2, src.cols - 2, CV_8UC1);for (int i = 1; i<src.rows - 1; i++) {for (int j = 1; j<src.cols - 1; j++) {//枚举在src中每个可以作为中心的像素位置_Tp center = src.at<_Tp>(i, j);unsigned char code = 0;//以i,j 为中心，从左上角顺时针对code进行编码code |= (src.at<_Tp>(i - 1, j - 1) > center) << 7;code |= (src.at<_Tp>(i - 1, j) > center) << 6;code |= (src.at<_Tp>(i - 1, j + 1) > center) << 5;code |= (src.at<_Tp>(i, j + 1) > center) << 4;code |= (src.at<_Tp>(i + 1, j + 1) > center) << 3;code |= (src.at<_Tp>(i + 1, j) > center) << 2;code |= (src.at<_Tp>(i + 1, j - 1) > center) << 1;code |= (src.at<_Tp>(i, j - 1) > center) << 0;dst.at<unsigned char>(i - 1, j - 1) = code;}}}//Extended LBP (Or Circular LBP)template <typename _Tp>void lbp::ELBP_(const Mat& src, Mat& dst, int radius, int neighbors) {//保证neighbors 的数量在1~31之间neighbors = max(min(neighbors, 31), 1); // set bounds...// Note: alternatively you can switch to the new OpenCV Mat_// type system to define an unsigned int matrix... I am probably// mistaken here, but I didn't see an unsigned int representation// in OpenCV's classic typesystem...//设置结果矩阵的行列数dst = Mat::zeros(src.rows - 2 * radius, src.cols - 2 * radius, CV_32SC1);//从x,y的取值来看，编码方向应该是从x正半轴逆时针旋转for (int n = 0; n < neighbors; n++) {// sample pointsfloat x = static_cast<float>(radius)* cos(2.0*M_PI*n / static_cast<float>(neighbors));float y = static_cast<float>(radius)* -sin(2.0*M_PI*n / static_cast<float>(neighbors));// relative indicesint 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));// fractional part 小数部分float ty = y - fy;float tx = x - fx;// set interpolation weightsfloat w1 = (1 - tx) * (1 - ty);float w2 = tx  * (1 - ty);float w3 = (1 - tx) *      ty;float w4 = tx  *      ty;// iterate through your datafor (int i = radius; i < src.rows - radius; i++) {for (int j = radius; j < src.cols - radius; j++) {//进行双线性插值float t = w1*src.at<_Tp>(i + fy, j + fx) + w2*src.at<_Tp>(i + fy, j + cx) + w3*src.at<_Tp>(i + cy, j + fx) + w4*src.at<_Tp>(i + cy, j + cx);// we are dealing with floating point precision, so add some little tolerancedst.at<unsigned int>(i - radius, j - radius) += ((t > src.at<_Tp>(i, j)) && (abs(t - src.at<_Tp>(i, j)) > std::numeric_limits<float>::epsilon())) << n;}}}//for int n neighbor}//这个没有看的很懂，因为感觉论文中的公式和下面和实现有点不对应template <typename _Tp>void lbp::VARLBP_(const Mat& src, Mat& dst, int radius, int neighbors) {max(min(neighbors, 31), 1); // set boundsdst = Mat::zeros(src.rows - 2 * radius, src.cols - 2 * radius, CV_32FC1); //! result// allocate some memory for temporary on-line variance calculationsMat _mean = Mat::zeros(src.rows, src.cols, CV_32FC1);Mat _delta = Mat::zeros(src.rows, src.cols, CV_32FC1);Mat _m2 = Mat::zeros(src.rows, src.cols, CV_32FC1);for (int n = 0; n < neighbors; n++) {// sample pointsfloat x = static_cast<float>(radius)* cos(2.0*M_PI*n / static_cast<float>(neighbors));float y = static_cast<float>(radius)* -sin(2.0*M_PI*n / static_cast<float>(neighbors));// relative indicesint 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));// fractional partfloat ty = y - fy;float tx = x - fx;// set interpolation weightsfloat w1 = (1 - tx) * (1 - ty);float w2 = tx  * (1 - ty);float w3 = (1 - tx) *      ty;float w4 = tx  *      ty;//和elbp比较相似，delta中保存着邻域点和中心点的差值，而mean中将差值按照一定的权值进行求和，先计算的邻域点权值越大//m2中值为 ： 本次的delta * （中心点像素值 - 本次mean值）// iterate through your datafor (int i = radius; i < src.rows - radius; i++) {for (int j = radius; j < src.cols - radius; j++) {float t = w1*src.at<_Tp>(i + fy, j + fx) + w2*src.at<_Tp>(i + fy, j + cx) + w3*src.at<_Tp>(i + cy, j + fx) + w4*src.at<_Tp>(i + cy, j + cx);_delta.at<float>(i, j) = t - _mean.at<float>(i, j);_mean.at<float>(i, j) = (_mean.at<float>(i, j) + (_delta.at<float>(i, j) / (1.0*(n + 1)))); // i am a bit paranoid_m2.at<float>(i, j) = _m2.at<float>(i, j) + _delta.at<float>(i, j) * (t - _mean.at<float>(i, j));}}}// calculate resultfor (int i = radius; i < src.rows - radius; i++) {for (int j = radius; j < src.cols - radius; j++) {dst.at<float>(i - radius, j - radius) = _m2.at<float>(i, j) / (1.0*(neighbors - 1));}}}// now the wrapper functionsvoid lbp::OLBP(const Mat& src, Mat& dst) {switch (src.type()) {case CV_8SC1: OLBP_<char>(src, dst); break;case CV_8UC1: OLBP_<unsigned char>(src, dst); break;case CV_16SC1: OLBP_<short>(src, dst); break;case CV_16UC1: OLBP_<unsigned short>(src, dst); break;case CV_32SC1: OLBP_<int>(src, dst); break;case CV_32FC1: OLBP_<float>(src, dst); break;case CV_64FC1: OLBP_<double>(src, dst); break;}}void lbp::ELBP(const Mat& src, Mat& dst, int radius, int neighbors) {switch (src.type()) {case CV_8SC1: ELBP_<char>(src, dst, radius, neighbors); break;case CV_8UC1: ELBP_<unsigned char>(src, dst, radius, neighbors); break;case CV_16SC1: ELBP_<short>(src, dst, radius, neighbors); break;case CV_16UC1: ELBP_<unsigned short>(src, dst, radius, neighbors); break;case CV_32SC1: ELBP_<int>(src, dst, radius, neighbors); break;case CV_32FC1: ELBP_<float>(src, dst, radius, neighbors); break;case CV_64FC1: ELBP_<double>(src, dst, radius, neighbors); break;}}void lbp::VARLBP(const Mat& src, Mat& dst, int radius, int neighbors) {switch (src.type()) {case CV_8SC1: VARLBP_<char>(src, dst, radius, neighbors); break;case CV_8UC1: VARLBP_<unsigned char>(src, dst, radius, neighbors); break;case CV_16SC1: VARLBP_<short>(src, dst, radius, neighbors); break;case CV_16UC1: VARLBP_<unsigned short>(src, dst, radius, neighbors); break;case CV_32SC1: VARLBP_<int>(src, dst, radius, neighbors); break;case CV_32FC1: VARLBP_<float>(src, dst, radius, neighbors); break;case CV_64FC1: VARLBP_<double>(src, dst, radius, neighbors); break;}}// now the Mat return functionsMat lbp::OLBP(const Mat& src) { Mat dst; OLBP(src, dst); return dst; }Mat lbp::ELBP(const Mat& src, int radius, int neighbors) { Mat dst; ELBP(src, dst, radius, neighbors); return dst; }Mat lbp::VARLBP(const Mat& src, int radius, int neighbors) { Mat dst; VARLBP(src, dst, radius, neighbors); return dst; }

2.2:lbpHistogram.cpp

#include "histogram.hpp"#include <vector>//numPatterns 为src中的模式数量template <typename _Tp>void lbp::histogram_(const Mat& src, Mat& hist, int numPatterns) {hist = Mat::zeros(1, numPatterns, CV_32SC1);for (int i = 0; i < src.rows; i++) {for (int j = 0; j < src.cols; j++) {int bin = src.at<_Tp>(i, j);hist.at<int>(0, bin) += 1;}}}template <typename _Tp>double lbp::chi_square_(const Mat& histogram0, const Mat& histogram1) {if (histogram0.type() != histogram1.type())CV_Error(CV_StsBadArg, "Histograms must be of equal type.");if (histogram0.rows != 1 || histogram0.rows != histogram1.rows || histogram0.cols != histogram1.cols)CV_Error(CV_StsBadArg, "Histograms must be of equal dimension.");double result = 0.0;for (int i = 0; i < histogram0.cols; i++) {double a = histogram0.at<_Tp>(0, i) - histogram1.at<_Tp>(0, i);double b = histogram0.at<_Tp>(0, i) + histogram1.at<_Tp>(0, i);if (abs(b) > numeric_limits<double>::epsilon()) {result += (a*a) / b;}}return result;}//创建空间直方图， 我觉得有bug//参数，hist为返回的空间直方图。numPatterns为每一个直方图的模式数量， window为在src上划分cell的大小//overlap为划分cell时，相邻cell重叠的大小//擦， 这个函数src的应该是对整幅图像处理过后的lbp特征矩阵，但这时候lbp特征矩阵和原始图像的大小已经不同，//这点需要注意void lbp::spatial_histogram(const Mat& src, Mat& hist, int numPatterns, const Size& window, int overlap) {int width = src.cols;int height = src.rows;vector<Mat> histograms;//个人觉得应该使用等号 原：x < width - window.width//y < height - window.heightfor (int x = 0; x < width - window.width; x += (window.width - overlap)) {for (int y = 0; y < height - window.height; y += (window.height - overlap)) {Mat cell = Mat(src, Rect(x, y, window.width, window.height));histograms.push_back(histogram(cell, numPatterns));}}hist.create(1, histograms.size()*numPatterns, CV_32SC1);// i know this is a bit lame now... feel free to make this a bit more efficient...//我倒是没啥办法能把它弄得more efficientfor (int histIdx = 0; histIdx < histograms.size(); histIdx++) {for (int valIdx = 0; valIdx < numPatterns; valIdx++) {int y = histIdx*numPatterns + valIdx;hist.at<int>(0, y) = histograms[histIdx].at<int>(valIdx);}}}// wrappersvoid lbp::histogram(const Mat& src, Mat& hist, int numPatterns) {switch (src.type()) {case CV_8SC1: histogram_<char>(src, hist, numPatterns); break;case CV_8UC1: histogram_<unsigned char>(src, hist, numPatterns); break;case CV_16SC1: histogram_<short>(src, hist, numPatterns); break;case CV_16UC1: histogram_<unsigned short>(src, hist, numPatterns); break;case CV_32SC1: histogram_<int>(src, hist, numPatterns); break;}}double lbp::chi_square(const Mat& histogram0, const Mat& histogram1) {switch (histogram0.type()) {case CV_8SC1: return chi_square_<char>(histogram0, histogram1); break;case CV_8UC1: return chi_square_<unsigned char>(histogram0, histogram1); break;case CV_16SC1: return chi_square_<short>(histogram0, histogram1); break;case CV_16UC1: return chi_square_<unsigned short>(histogram0, histogram1); break;case CV_32SC1: return chi_square_<int>(histogram0, histogram1); break;}}void lbp::spatial_histogram(const Mat& src, Mat& dst, int numPatterns, int gridx, int gridy, int overlap) {int width = static_cast<int>(floor(src.cols / gridx));int height = static_cast<int>(floor(src.rows / gridy));spatial_histogram(src, dst, numPatterns, Size_<int>(width, height), overlap);}// Mat return type functionsMat lbp::histogram(const Mat& src, int numPatterns) {Mat hist;histogram(src, hist, numPatterns);return hist;}Mat lbp::spatial_histogram(const Mat& src, int numPatterns, const Size& window, int overlap) {Mat hist;spatial_histogram(src, hist, numPatterns, window, overlap);return hist;}Mat lbp::spatial_histogram(const Mat& src, int numPatterns, int gridx, int gridy, int overlap) {Mat hist;spatial_histogram(src, hist, numPatterns, gridx, gridy);return hist;}

参考文献：

【1】. Ojala T, Pietikainen M, Maenpaa T. Multiresolution gray-scale and rotation invariant texture classification with local binary patterns[J]. Pattern Analysis and Machine Intelligence, IEEE Transactions on, 2002, 24(7): 971-987.