老外的原文：《Multiresolution gray-scale and rotation invariant texture classification with local binary patterns》 Timo Ojala, Matti PietikaÈ inen,2002
本文将对这篇文章做部分翻译，最后将分别给出下列代码实现
。灰度不变性LBP( gray scale invariant)
。旋转不变性LBP(rotation invariant)
。旋转不变等价LBP(rotation & uniform invariant)
1.INTRODUCTION
二维纹理分析有许多潜在的应用，比如工业表面检查、遥感、生物图片分析，但是只有很少的例子成功发现了纹理的存在。主要的问题是，在现实世界中，由于方向、尺度及其他视觉外观的变化导致纹理常常不是均一的。灰度尺度不变性通常变得非常重要，因为不均匀的光照或图像内大幅变化。另外，已经提出的纹理测量算法计算量太大【比如灰度共生矩阵，我加的】。
绝大多数纹理分类方法都显式或隐式的假设，待分类的未知样本与训练样本在在空间大小、方向和灰度尺度属性上都一致。然而，现实世界的纹理受变化的光照条件影响，可以存在任意空间分辨率和旋转方向。
然后，罗列了一堆前人的工作......
在老外这篇文章中，老外提出了一种计算量小并且可以有效检测大量旋转和灰度尺度不同的纹理。他们提出一种基于Local Binary Patterns【LBP特征】 的灰度尺度不变性和旋转不变性的纹理算子。从局部邻域的圆形区域内像素灰度的联合分布开始，他们推导了一个算子，对任何单调的灰度尺度变化都能显示不变性。旋转不变性是通过灰度尺度不变算子配合一系列旋转不变模式实现的。
【我这算直译的，不通也是常理，下面开始重点了】
2.GRAY SCALE AND ROTATION INVARIANT LOCAL BINARY PATTERNS
定义纹理 T 为一副黑白纹理图像局部邻域内图像像素灰度的联合分布：
其中，gc 对应局部邻域中心像素的灰度值，gp (p = 1,2,...,P-1)对应半径为 R 的圆形对称邻域内 P等分的像素灰度。如果 gc 坐标为 (0,0)，那么 gp 的坐标为…（- R sin…（2πp/P），R cos… （2πp/P））†。如下图所示，为不同采样点 P 和半径 R 的圆形对称邻域。其中对于邻域内那些没有直接落在像素【方格】*的点的灰度值将通过二线性插值完成。
2.1 Achieving Gray-Scale Invariance
首先，我们将邻域内中心像素点 gc 的灰度值与圆形对称区域内其他点 gp 的灰度值相减：
接下来，我们假设差分 gp - gc 和 gc 是相互独立的，于是就有：
在实际中，这种直接独立是不能被保证的，因此上式仅是联合分布的一种近似。由于分布 t ( gc )描述了图像整体光照，因此没有给纹理分析提供有用的信息。这样，纹理特征可以由以下联合差分分布表示：
这是一个具有高度判别性的纹理算子。它在 P 维直方图中，记录了每个像素点的邻域中模式变化的发生次数。对于恒定的区域，所有方向的差分均为 0。对于缓慢的斜边，该算子记录了梯度方向最大的差分，沿着边的方向差分为0 。对于一个斑点，所有方向的差分都很大。
带符号差分 gp - gc 不受平均光源变化影响。因此，联合差分分布对于灰度尺度变化具有不变性。我们通过仅使用差分的符号，而不是差分本身的值来实现灰度尺度不变性：
通过对每个符号 s( gp - gc ) 赋予因子2^P，我们将上述公式变成唯一的 LBP值，用来描述局部图像纹理的空间结构：
当P=8,R=1时，这里的LBP值和经典的LBP值很相似，但差别在于，这边采样区域为圆形，且像素值必须通过二线性插值获得。
2.2 Achieving Rotation Invariance
上述LBP算子将产生 2^P 个不同输出，对应由邻域内P个像素点形成的2 ^ P个不同二进制模式。当图像被旋转，【以下为意译】，灰度值gp 和g0的相对位置将发生变化，而g0通常取中心 gc 的正右 ，坐标(0，R) ，这会导致不同的LBP值。但任何角度的旋转不影响，圆形对称邻域内 0 和 1 的相对位置关系，为了移除旋转获得唯一的LBP值，我们定义：
其中，ROR（x，i）执行沿时钟方向将P位数 X 移动 i 次。对于图像像素而言，就是将邻域集合按照时钟方向旋转很多次，直到当前旋转下构成的LBP值最小。【即文中所说的，最重要位的最大值为0，即gp-1为0】
如下图所示，对于8个采样点，将有36种唯一的旋转不变二值模式：【黑点为0，白点为1】
2.3 Improved Rotation Invariance with ªUniformº Patterns
对于8个采样点，灰度不变性LBP将产生256中输出，旋转不变性LBP将产生36个输出，而基于unifrom的旋转不变LBP将只有9中输出。【uniform形式有58中输出】
【上面讲的种类对于后面编程实现，计算映射mapping很重要】
首先介绍什么是uniform，它是指均匀环形结构内包含非常少的空间转换。我们定义U(pattern),用来记录空间转换的数量，即0-1变化的次数。比如，上图第一行中00000000和11111111，它们的U值均为0，而其余种类的U均为2.类似的，上述其他27中模式的 U 值至少为4。我们将 U 值为2的模式称为 uniform
那么：
对于8个采样点，共有58种唯一的 uniform类型的LBP值输出【此处没考虑旋转不变性，58种输出是LBP的2进制值，而非U值，U值都是2】，如下图所示：
关于uniform输出的种类，并不在该论文中，属于题外话，因为有人直接拿uniform形式去提取LBP特征。详情见《Rotation Invariant Image Description with Local Binary Pattern Histogram Fourier Features》一文。
接下来继续介绍基于uniform的旋转不变LBP，我们定义如下公式：
其中，P为采样点个数。
3.具体实现
下面代码用C实现，读入一副灰度图像，采样点个数，采样半径，对每个像素计算LBP特征并输出，LBP图像及其直方图。
对于旋转不变性及uniform旋转不变性，我没有对每个像素都进行上述公式的操作，由于是2进制的圆形循环，可以提前做个mapping映射关系，加快程序执行速度。
采样点数量均为8，采样半径为10
输入图像：
3.1 灰度不变性LBP：
void gray_invariant_lbp(IplImage *src,int height,int width,int num_sp,MyPoint *spoint)
{
IplImage *target,*hist;
int i,j,k,box_x,box_y,orign_x,orign_y,dx,dy,tx,ty,fy,fx,cy,cx,v;
double min_x,max_x,min_y,max_y,w1,w2,w3,w4,N,x,y;
int *result;
float dishu;
dishu = 2.0;
max_x=0;max_y=0;min_x=0;min_y=0;
for (k=0;k<num_sp;k++)
{
if (max_x<spoint[k].x)
{
max_x=spoint[k].x;
}
if (max_y<spoint[k].y)
{
max_y=spoint[k].y;
}
if (min_x>spoint[k].x)
{
min_x=spoint[k].x;
}
if (min_y>spoint[k].y)
{
min_y=spoint[k].y;
}
}
//计算模版大小
box_x = ceil(MAX(max_x,0)) - floor(MIN(min_x,0)) + 1;
box_y = ceil(MAX(max_y,0)) - floor(MIN(min_y,0)) + 1;
if (width<box_x||height<box_y)
{
printf("Too small input image. Should be at least (2*radius+1) x (2*radius+1)");
return;
}
//计算可滤波图像大小,opencv图像数组下标从0开始
orign_x = 0 - floor(MIN(min_x,0));//起点
orign_y = 0 - floor(MIN(min_x,0));
dx = width - box_x+1;//终点
dy = height - box_y+1;
int cols = pow(dishu,(float)num_sp);
hist = cvCreateImage(cvSize(300,200),IPL_DEPTH_8U,3);//直方图图像
target = cvCreateImage(cvSize(dx,dy),IPL_DEPTH_8U,1);
result = (int *)malloc(sizeof(int)*dx*dy);
double *val_hist = (double *)malloc(sizeof(double)*cols); //直方图数组
memset(result,0,sizeof(int)*dx*dy);
CvRect roi =cvRect(orign_x, orign_y, dx, dy);
cvSetImageROI(src, roi);
cvCopy(src, target);
cvResetImageROI(src);
cvSaveImage("haha.jpg",target);
for ( k = 0; k<num_sp;k++)
{
x = spoint[k].x+orign_x;
y = spoint[k].y+orign_y;
//二线性插值图像
fy = floor(y); //向下取整
fx = floor(x);
cy = ceil(y); //向上取整
cx = ceil(x);
ty = y - fy;
tx = x - fx;
w1 = (1 - tx) * (1 - ty);
w2 = tx * (1 - ty);
w3 = (1 - tx) * ty ;
w4 = tx * ty ;
v = pow(dishu,(float)k);
for (i = 0;i<dy;i++)
{
for (j = 0;j<dx;j++)
{
//灰度插值图像像素
N = w1 * (double)(unsigned char)src->imageData[(i+fy)*src->width+j+fx]+
w2 * (double)(unsigned char)src->imageData[(i+fy)*src->width+j+cx]+
w3 * (double)(unsigned char)src->imageData[(i+cy)*src->width+j+fx]+
w4 * (double)(unsigned char)src->imageData[(i+cy)*src->width+j+cx];
if( N >= (double)(unsigned char)target->imageData[i*dx+j])
{
result[i*dx+j] = result[i*dx+j] + v * 1;
}else{
result[i*dx+j] = result[i*dx+j] + v * 0;
}
}
}
}
//显示图像
if (num_sp<=8)
{
//只有采样数小于8，则编码范围0-255，才能显示图像
for (i = 0;i<dy;i++)
{
for (j = 0;j<dx;j++)
{
target->imageData[i*dx+j] = (unsigned char)result[i*dx+j];
//printf("%d\n",(unsigned char)target->imageData[i*width+j]);
}
}
cvSaveImage("result.jpg",target);
}
//显示直方图
for (i=0;i<cols;i++)
{
val_hist[i]=0.0;
}
for (i=0; i<dy*dx;i++)
{
val_hist[result[i]]++;
}
double temp_max=0.0;
for (i=0;i<cols;i++) //求直方图最大值，为了归一化
{
//printf("%f\n",val_hist[i]);
if (temp_max<val_hist[i])
{
temp_max=val_hist[i];
}
}
//画直方图
CvPoint p1,p2;
double bin_width=(double)hist->width/cols;
double bin_unith=(double)hist->height/temp_max;
for (i=0;i<cols;i++)
{
p1.x=i*bin_width;p1.y=hist->height;
p2.x=(i+1)*bin_width;p2.y=hist->height-val_hist[i]*bin_unith;
cvRectangle(hist,p1,p2,cvScalar(0,255),-1,8,0);
}
cvSaveImage("hist.jpg",hist);
}
实验结果：
3.2 旋转不变性：
void rotation_invariant_mapping(int range,int num_sp,int *Mapping)
{
int newMax,rm,r;
int *tmpMap;
newMax = 0;
tmpMap = (int *)malloc(sizeof(int)*range);
memset(tmpMap,-1,sizeof(int)*range);
for (int i = 0 ; i < range ; i++)
{
rm = i;
r = i;
for (int j = 0 ; j < num_sp -1 ;j++)
{
//将r向左循环移动一位,当r超过num_sp位时，舍弃
r = r << 1;
if (r > range -1)
{
r = r - (range -1);
}
//printf("%d,%d\n",r,rm);
if (r < rm)
{
rm = r;
}
}
if (tmpMap[rm] < 0)
{
tmpMap[rm] = newMax;
newMax++;
}
Mapping[i] = tmpMap[rm];
}
free(tmpMap);
}
void rotation_invariant_lbp(IplImage *src,int height,int width,int num_sp,MyPoint *spoint,int *Mapping)
{
IplImage *target,*hist;
int i,j,k,box_x,box_y,orign_x,orign_y,dx,dy,tx,ty,fy,fx,cy,cx,v;
double min_x,max_x,min_y,max_y,w1,w2,w3,w4,N,x,y;
int *result;
float dishu;
dishu = 2.0;
max_x=0;max_y=0;min_x=0;min_y=0;
for (k=0;k<num_sp;k++)
{
if (max_x<spoint[k].x)
{
max_x=spoint[k].x;
}
if (max_y<spoint[k].y)
{
max_y=spoint[k].y;
}
if (min_x>spoint[k].x)
{
min_x=spoint[k].x;
}
if (min_y>spoint[k].y)
{
min_y=spoint[k].y;
}
}
//计算模版大小
box_x = ceil(MAX(max_x,0)) - floor(MIN(min_x,0)) + 1;
box_y = ceil(MAX(max_y,0)) - floor(MIN(min_y,0)) + 1;
if (width<box_x||height<box_y)
{
printf("Too small input image. Should be at least (2*radius+1) x (2*radius+1)");
return;
}
//计算可滤波图像大小,opencv图像数组下标从0开始
orign_x = 0 - floor(MIN(min_x,0));//起点
orign_y = 0 - floor(MIN(min_x,0));
dx = width - box_x+1;//终点
dy = height - box_y+1;
target = cvCreateImage(cvSize(dx,dy),IPL_DEPTH_8U,1);
result = (int *)malloc(sizeof(int)*dx*dy);
memset(result,0,sizeof(int)*dx*dy);
CvRect roi =cvRect(orign_x, orign_y, dx, dy);
cvSetImageROI(src, roi);
cvCopy(src, target);
cvResetImageROI(src);
cvSaveImage("haha.jpg",target);
for ( k = 0; k<num_sp;k++)
{
x = spoint[k].x+orign_x;
y = spoint[k].y+orign_y;
//二线性插值图像
fy = floor(y); //向下取整
fx = floor(x);
cy = ceil(y); //向上取整
cx = ceil(x);
ty = y - fy;
tx = x - fx;
w1 = (1 - tx) * (1 - ty);
w2 = tx * (1 - ty);
w3 = (1 - tx) * ty ;
w4 = tx * ty ;
v = pow(dishu,(float)k);
for (i = 0;i<dy;i++)
{
for (j = 0;j<dx;j++)
{
//灰度插值图像像素
N = w1 * (double)(unsigned char)src->imageData[(i+fy)*src->width+j+fx]+
w2 * (double)(unsigned char)src->imageData[(i+fy)*src->width+j+cx]+
w3 * (double)(unsigned char)src->imageData[(i+cy)*src->width+j+fx]+
w4 * (double)(unsigned char)src->imageData[(i+cy)*src->width+j+cx];
if( N >= (double)(unsigned char)target->imageData[i*dx+j])
{
result[i*dx+j] = result[i*dx+j] + v * 1;
}else{
result[i*dx+j] = result[i*dx+j] + v * 0;
}
}
}
}
//将result的值映射为mapping的值
for(i = 0; i < dy ;i++)
{
for (j = 0; j < dx ;j ++)
{
result[i*dx+j] = Mapping[result[i*dx+j]];
}
}
//显示图像
int cols = 0;//直方图的横坐标，也是result数组的元素种类
int mapping_size = pow(dishu,(float)num_sp);
for (i = 0;i < mapping_size; i++ )
{
if (cols < Mapping[i])
{
cols = Mapping[i];
}
}
if (cols < 255)
{
//只有采样数小于8，则编码范围0-255，才能显示图像
for (i = 0;i<dy;i++)
{
for (j = 0;j<dx;j++)
{
target->imageData[i*dx+j] = (unsigned char)result[i*dx+j];
//printf("%d\n",(unsigned char)target->imageData[i*width+j]);
}
}
cvSaveImage("result.jpg",target);
}
//计算直方图
hist = cvCreateImage(cvSize(300,200),IPL_DEPTH_8U,3);//直方图图像
double *val_hist = (double *)malloc(sizeof(double)*cols); //直方图数组
for (i=0;i<cols;i++)
{
val_hist[i]=0.0;
}
for (i=0; i<dy*dx;i++)
{
val_hist[result[i]]++;
}
double temp_max=0.0;
for (i=0;i<cols;i++) //求直方图最大值，为了归一化
{
//printf("%f\n",val_hist[i]);
if (temp_max<val_hist[i])
{
temp_max=val_hist[i];
}
}
//画直方图
CvPoint p1,p2;
double bin_width=(double)hist->width/cols;
double bin_unith=(double)hist->height/temp_max;
for (i=0;i<cols;i++)
{
p1.x=i*bin_width;p1.y=hist->height;
p2.x=(i+1)*bin_width;p2.y=hist->height-val_hist[i]*bin_unith;
cvRectangle(hist,p1,p2,cvScalar(0,255),-1,8,0);
}
cvSaveImage("hist.jpg",hist);
}
实验结果：
3.3 uniform 旋转不变性，给出完整代码：
#include "stdafx.h"
#include "cv.h"
#include "highgui.h"
#define PI 3.1415926
#define MAX(x,y) (x)>(y)?(x):(y)
#define MIN(x,y) (x)<(y)?(x):(y)
typedef struct MyPoint
{
double x;
double y;
}MyPoint;
void calc_position(int radius,int num_sp,MyPoint *spoint)
{
double theta;
theta = 2*PI/num_sp;
for (int i = 0; i < num_sp; i++)
{
spoint[i].y = -radius * sin(i * theta);
spoint[i].x = radius * cos(i * theta);
}
}
int calc_sum(int r)
{
int res_sum;
res_sum = 0;
while (r)
{
res_sum = res_sum + r % 2;
r /= 2;
}
return res_sum;
}
void rotation_uniform_invariant_mapping(int range,int num_sp,int *Mapping)
{
int numt,i,j,tem_xor;
numt = 0;
tem_xor = 0;
for (i = 0; i< range; i++)
{
j = i << 1;
if (j > range -1)
{
j = j - (range -1);
}
tem_xor = i ^ j; // 异或
numt = calc_sum(tem_xor);//计算异或结果中1的个数，即跳变个数
if (numt <= 2)
{
Mapping[i] = calc_sum(i);
}else{
Mapping[i] = num_sp+1;
}
}
}
void rotation_uniform_invariant_lbp(IplImage *src,int height,int width,int num_sp,MyPoint *spoint,int *Mapping)
{
IplImage *target,*hist;
int i,j,k,box_x,box_y,orign_x,orign_y,dx,dy,tx,ty,fy,fx,cy,cx,v;
double min_x,max_x,min_y,max_y,w1,w2,w3,w4,N,x,y;
int *result;
float dishu;
dishu = 2.0;
max_x=0;max_y=0;min_x=0;min_y=0;
for (k=0;k<num_sp;k++)
{
if (max_x<spoint[k].x)
{
max_x=spoint[k].x;
}
if (max_y<spoint[k].y)
{
max_y=spoint[k].y;
}
if (min_x>spoint[k].x)
{
min_x=spoint[k].x;
}
if (min_y>spoint[k].y)
{
min_y=spoint[k].y;
}
}
//计算模版大小
box_x = ceil(MAX(max_x,0)) - floor(MIN(min_x,0)) + 1;
box_y = ceil(MAX(max_y,0)) - floor(MIN(min_y,0)) + 1;
if (width<box_x||height<box_y)
{
printf("Too small input image. Should be at least (2*radius+1) x (2*radius+1)");
return;
}
//计算可滤波图像大小,opencv图像数组下标从0开始
orign_x = 0 - floor(MIN(min_x,0));//起点
orign_y = 0 - floor(MIN(min_x,0));
dx = width - box_x+1;//终点
dy = height - box_y+1;
target = cvCreateImage(cvSize(dx,dy),IPL_DEPTH_8U,1);
result = (int *)malloc(sizeof(int)*dx*dy);
memset(result,0,sizeof(int)*dx*dy);
CvRect roi =cvRect(orign_x, orign_y, dx, dy);
cvSetImageROI(src, roi);
cvCopy(src, target);
cvResetImageROI(src);
cvSaveImage("haha.jpg",target);
for ( k = 0; k<num_sp;k++)
{
x = spoint[k].x+orign_x;
y = spoint[k].y+orign_y;
//二线性插值图像
fy = floor(y); //向下取整
fx = floor(x);
cy = ceil(y); //向上取整
cx = ceil(x);
ty = y - fy;
tx = x - fx;
w1 = (1 - tx) * (1 - ty);
w2 = tx * (1 - ty);
w3 = (1 - tx) * ty ;
w4 = tx * ty ;
v = pow(dishu,(float)k);
for (i = 0;i<dy;i++)
{
for (j = 0;j<dx;j++)
{
//灰度插值图像像素
N = w1 * (double)(unsigned char)src->imageData[(i+fy)*src->width+j+fx]+
w2 * (double)(unsigned char)src->imageData[(i+fy)*src->width+j+cx]+
w3 * (double)(unsigned char)src->imageData[(i+cy)*src->width+j+fx]+
w4 * (double)(unsigned char)src->imageData[(i+cy)*src->width+j+cx];
if( N >= (double)(unsigned char)target->imageData[i*dx+j])
{
result[i*dx+j] = result[i*dx+j] + v * 1;
}else{
result[i*dx+j] = result[i*dx+j] + v * 0;
}
}
}
}
//将result的值映射为mapping的值
for(i = 0; i < dy ;i++)
{
for (j = 0; j < dx ;j ++)
{
result[i*dx+j] = Mapping[result[i*dx+j]];
}
}
//显示图像
int cols = 0;//直方图的横坐标，也是result数组的元素种类
int mapping_size = pow(dishu,(float)num_sp);
for (i = 0;i < mapping_size; i++ )
{
if (cols < Mapping[i])
{
cols = Mapping[i];
}
}
if (cols < 255)
{
//只有采样数小于8，则编码范围0-255，才能显示图像
for (i = 0;i<dy;i++)
{
for (j = 0;j<dx;j++)
{
target->imageData[i*dx+j] = (unsigned char)result[i*dx+j];
//printf("%d\n",(unsigned char)target->imageData[i*width+j]);
}
}
cvSaveImage("result.jpg",target);
}
//计算直方图
hist = cvCreateImage(cvSize(300,200),IPL_DEPTH_8U,3);//直方图图像
double *val_hist = (double *)malloc(sizeof(double)*cols); //直方图数组
for (i=0;i<cols;i++)
{
val_hist[i]=0.0;
}
for (i=0; i<dy*dx;i++)
{
val_hist[result[i]]++;
}
double temp_max=0.0;
for (i=0;i<cols;i++) //求直方图最大值，为了归一化
{
//printf("%f\n",val_hist[i]);
if (temp_max<val_hist[i])
{
temp_max=val_hist[i];
}
}
//画直方图
CvPoint p1,p2;
double bin_width=(double)hist->width/cols;
double bin_unith=(double)hist->height/temp_max;
for (i=0;i<cols;i++)
{
p1.x=i*bin_width;p1.y=hist->height;
p2.x=(i+1)*bin_width;p2.y=hist->height-val_hist[i]*bin_unith;
cvRectangle(hist,p1,p2,cvScalar(0,255),-1,8,0);
}
cvSaveImage("hist.jpg",hist);
}
int _tmain(int argc, _TCHAR* argv[])
{
IplImage *src,*grey,*result;
int samples,radius,range,*mapping;
MyPoint *spoint;
float Mi;
samples = 8;
radius = 10;
Mi = 2.0;
range = pow(Mi,samples);
src = cvLoadImage("test2.jpg");
grey = cvCreateImage(cvSize(src->width,src->height),IPL_DEPTH_8U,1);
cvCvtColor(src,grey,CV_BGR2GRAY);
mapping = (int *)malloc(sizeof(int)*range);
memset(mapping,0,sizeof(int)*range);
//计算采样点相对坐标
spoint = (MyPoint *)malloc(sizeof(MyPoint)*samples);
calc_position(radius,samples,spoint);
//计算灰度不变性LBP特征，写回浮点数图像矩阵中
//gray_invariant_lbp(grey,src->height,src->width,samples,spoint);
//计算旋转不变形LBP特征
//rotation_invariant_mapping(range,samples,mapping);
//rotation_invariant_lbp(grey,src->height,src->width,samples,spoint,mapping);
//计算旋转不变等价LBP特征
rotation_uniform_invariant_mapping(range,samples,mapping);
rotation_uniform_invariant_lbp(grey,src->height,src->width,samples,spoint,mapping);
return 0;
}
实验结果：
越到后面越黑，因为输出的种类越来越少。