python实现LBP纹理提取

• 什么是LBP纹理特征？
  LBP（Local Binary Patterns，局部二值模式）是提取局部特征作为判别依据的，一种有效的纹理描述算子，度量和提取图像局部的纹理信息。它具有旋转不变性和灰度不变性等显著的优点，对光照具有不变性。由T. Ojala, M.Pietikäinen, 和D. Harwood 在1994年提出用于纹理特征提取。
  LBP特征基本描述：
  LBP的基本思想是定义于像素的8邻域中（3x3的窗口）, 以中心像素的灰度值为阈值, 将周围8 个像素的值与其比较, 如果周围的像素值小于中心像素的灰度值, 该像素位置就被标记为0, 否则标记为1。这样3x3的邻域内的八个点经过比较，顺时针能够产生8位二进制数（通常转换为十进制，即LBP码，共256种），每个像素得到一个二进制组合, 即得到该窗口中心像素点的LBP值，并用这个值来反映该区域的纹理信息，如下图。
  
  LBP算法的核心思想，是以某个像素点为中心，与其邻域像素点共同计算。所以光照变化很难改变中心像素点的灰度值与周围8个像素的大小关系。因为光照是一种区域性质的变化，而不是单像素性质的变化。
  LBP特征改进描述：
  1）圆形LBP算子：
  基本的 LBP算子的最大缺陷在于它只覆盖了一个$3\ast 3$固定范围内的小区域，这显然不能满足不同尺寸纹理的需要。为了适应不同尺寸的纹理特征，达到灰度和旋转不变性的要求。Ojala将 3×3邻域扩展到任意邻域，并用圆形邻域代替了正方形邻域。改进后的 LBP 算子允许在半径为 R的圆形邻域内有任意多个像素点。从而得到了半径为R的圆形区域内含有P个采样点的LBP算子，如下图。（2）2）LBP旋转不变模式算子：
  为保证图像旋转时LBP不变，提出了LBP旋转不变算法。即不断旋转圆形邻域得到一系列初始定义的 LBP值，取其最小值作为该邻域的 LBP 值。这样就出现了多个位置模式对应同一个LBP值的情况，但这其实保证了图像旋转以后还是同一个值。如下图，8种模式都对应的一个LBP(00001111)，8种模式可以看做一个图旋转个不同的角度（黑色0，白色1）。
  
  3）LBP等价算子：
  由上可知LBP会产生很多不同的模式，比如半径5,20个采样点，就有可能产生$2^{20}$种二进制模式，这对图像整体纹理描述很不方便。即使在纹理分类或人脸识别时采用的是LBP模式的统计直方图来进行的，但数据量过大将不合适，且直方图表示过于稀疏。所以需要对LBP模式进行降维，使得利用更少的数据就可以表达很好的图像纹理信息。
  Ojala提出了（Uniform Pattern）进行降维，利用LBP算子中0,1改变的次数来进行归类，定义从0到1或从1到0最多有两次突变时，该LBP对应的称为一个等价模式，而等价模式以为的都归为混合模式。这样等价模式的总数为2+2(p-1)+(p-1)(p-2)/2种，其中p为采样点数。0次突变2种，1次突变2(p-1)种，两次突变(p-1)(p-2)/2种。
• LBP算子怎么用python实现？
  采用$skimage$包中的local_binary_pattern
  skimage.feature.local_binary_pattern(image, P, R, method='default')
  灰度和旋转不变的LBP(局部二进制模式)。
  参数：
  image：(N, M) array
  灰度图像
  P：int
  圆对称邻居设定点的数目(角空间的量化)（采样点的数量）。
  R：float
  圆的半径(算子的空间分辨率)（圆形算子的半径）。
  method：{‘default’, ‘ror’, ‘uniform’, ‘var’}
  确定算子计算的模式。
  ‘default’: 原始的局部二值模式是灰度的，但不是旋转不变（基本算子）。
  ‘ror’: 扩展的模式（旋转不变模式算子）实现是灰度和旋转不变。
  ‘uniform’: 改进的旋转不变性均匀模式和对具有灰度和旋转不变性的角空间进行精细量化（降维算子）。
  ‘nri_uniform’: 非旋转不变的均匀模式变量，只有灰度不变（非旋转+降维）
  ‘var’: 旋转不变性方差是对局部图像纹理对比度的度量，它具有旋转不变性而不具有灰度不变性。
  返回：
  output：(N, M) array
  LBP 图像.
  LBP算子对于图片的实际物理意义理解：
  
  上图可以理解为：黑色(或白色)表示比中心像素更弱(或更强)的像素。当周围的像素都是黑色或白色，那么图像区域是平坦的(即无特征)。连续的黑色或白色像素组被认为是统一的图案，可以解释为角或边。如果像素在黑白像素之间来回切换，则该模式被认为是不均匀的。
• 应用实例：

from skimage.feature import local_binary_pattern
from skimage import data
from skimage.color import label2rgb
import matplotlib.pyplot as plt
import numpy as np
radius = 3
n_points = 8 * radius
image = data.brick()
plt.imshow(image,cmap='gray')
METHOD = 'uniform'#采用等效算子
lbp = local_binary_pattern(image, n_points, radius, METHOD)
plt.imshow(lbp,cmap='gray')
plt.colorbar()
def hist(ax, lbp):#绘制LBP的直方图
    n_bins = int(lbp.max() + 1)#直方图横坐标
    return ax.hist(lbp.ravel(), density=True, bins=n_bins, range=(0, n_bins),
                   facecolor='0.5')

w = width = radius - 1
edge_labels = range(n_points // 2 - w, n_points // 2 + w + 1)#边点对应的LBP范围
flat_labels = list(range(0, w + 1)) + list(range(n_points - w, n_points + 2))#平面点对应的LBP范围
i_14 = n_points // 4            #直方图的四分之一
i_34 = 3 * (n_points // 4) #直方图的四分之三
corner_labels = (list(range(i_14 - w, i_14 + w + 1)) +
                 list(range(i_34 - w, i_34 + w + 1)))#角点对应的LBP范围
label_sets = (edge_labels, flat_labels, corner_labels)#三种点作为标签集
def overlay_labels(image, lbp, labels):
    mask = np.logical_or.reduce([lbp == each for each in labels])
    return label2rgb(mask, image=image, bg_label=0, alpha=0.5)
    #返回一个RGB图像，其中彩色编码标签被涂在图像上
def highlight_bars(bars, indexes):
    for i in indexes:
        bars[i].set_facecolor('r')
        #在直方图里将直方图的条状高亮为红色      
fig, (ax_img, ax_hist) = plt.subplots(nrows=2, ncols=3, figsize=(9, 6))
plt.gray()
titles = ('edge', 'flat', 'corner')
for ax, labels in zip(ax_img, label_sets):
    ax.imshow(overlay_labels(image, lbp, labels))#在原图上绘制标记点
for ax, labels, name in zip(ax_hist, label_sets, titles):
    counts, _, bars = hist(ax, lbp)
    highlight_bars(bars, labels)
    ax.set_ylim(top=np.max(counts[:-1]))
    ax.set_xlim(right=n_points + 2)
    ax.set_title(name)
ax_hist[0].set_ylabel('Percentage')#设置y标签为百分数
for ax in ax_img:
    ax.axis('off')

LBP算子经常与直方图联合使用:

radius = 2
n_points = 8 * radius


def kullback_leibler_divergence(p, q):
    p = np.asarray(p)
    q = np.asarray(q)
    filt = np.logical_and(p != 0, q != 0)
    return np.sum(p[filt] * np.log2(p[filt] / q[filt]))


def match(refs, img):
    best_score = 10
    best_name = None
    lbp = local_binary_pattern(img, n_points, radius, METHOD)
    n_bins = int(lbp.max() + 1)
    hist, _ = np.histogram(lbp, density=True, bins=n_bins, range=(0, n_bins))
    for name, ref in refs.items():
        ref_hist, _ = np.histogram(ref, density=True, bins=n_bins,
                                   range=(0, n_bins))
        score = kullback_leibler_divergence(hist, ref_hist)
        if score < best_score:
            best_score = score
            best_name = name
    return best_name


brick = data.brick()
grass = data.grass()
gravel = data.gravel()

refs = {
    'brick': local_binary_pattern(brick, n_points, radius, METHOD),
    'grass': local_binary_pattern(grass, n_points, radius, METHOD),
    'gravel': local_binary_pattern(gravel, n_points, radius, METHOD)
}

# classify rotated textures
print('Rotated images matched against references using LBP:')
print('original: brick, rotated: 30deg, match result: ',
      match(refs, rotate(brick, angle=30, resize=False)))
print('original: brick, rotated: 70deg, match result: ',
      match(refs, rotate(brick, angle=70, resize=False)))
print('original: grass, rotated: 145deg, match result: ',
      match(refs, rotate(grass, angle=145, resize=False)))

# plot histograms of LBP of textures
fig, ((ax1, ax2, ax3), (ax4, ax5, ax6)) = plt.subplots(nrows=2, ncols=3,
                                                       figsize=(9, 6))
plt.gray()

ax1.imshow(brick)
ax1.axis('off')
hist(ax4, refs['brick'])
ax4.set_ylabel('Percentage')

ax2.imshow(grass)
ax2.axis('off')
hist(ax5, refs['grass'])
ax5.set_xlabel('Uniform LBP values')

ax3.imshow(gravel)
ax3.axis('off')
hist(ax6, refs['gravel'])

plt.show()

程序1结果：

程序二结果：