介绍

机器学习中特征提取的共同目标是将原始数据表示为一组简化的特征，以便更好地描述其主要特征和属性。通过这种方式，我们可以降低机器学习中原始输入的维数，并使用新特征作为输入来训练模式识别和分类技术。

虽然我们可以从图片中提取出一些特征，但是局部二进制模式(LBP)在理论上是一种简单而有效的灰度和旋转不变纹理分类方法。它们之所以有效，是因为最常见的模式对应于原始的微特征，如边缘、角落、斑点、平坦区域。

Ojala等研究表明，均匀图案的离散发生直方图是一个非常强大的纹理特征。图像纹理是一种二维现象，具有两种特性:(1)空间结构(图案)和(2)对比度。

方法

Circularly Symmetric Neighbor Set

对于给定的像素gc，Circularly Symmetric Neighbor Set是由在半径为R的圆上围绕中心点的坐标点(i, j)和一些元素P定义的。

纹理

我们将纹理T定义为灰度图像中的像素集合

其中gp对应于p局部邻域的灰度值。

插值

当邻域不在像素的中心时，应该通过插值计算该邻域灰度值。因此，我们需要定义一个给定坐标的函数，返回插值的灰度值。

实现灰度不变性

考虑到可能的信息丢失，可以将纹理转换为joint difference。为了计算它，我们将中心像素的灰度值减去所有的neighbor set。联合差分分布是一种高度区分的纹理算子。它记录了在p维直方图中每个像素附近出现的各种图案。

其中gp是p相邻的灰色值。这种分布对于灰度变化是不变的。

局部二进制模式

根据定义，LBP_ {P，R}运算符对于灰度的任何单调变换是不变的。只要灰度值的顺序保持不变，LBP_ {P，R}运算符的输出保持不变。

其中

Uniform 局部二进制模式

Ojala提到，在他们的实践经验中，LBP不是一个好的判别器。他们建议只选择一组局部二进制模式，使空间过渡的数量（按位0/1变化）不超过2.例如，模式'1111'有0个空间过渡，模式'1100'有1个空间过渡和模式'1101'具有2个空间过渡。对于每个统一模式，关联唯一索引。

现在，我们可以计算中心像素的局部二进制模式。下一步是计算所有像素的局部二进制模式。

提示：为简单起见，我不考虑所选索引为负的情况（即img_gray [-1] [0]返回第一列的最后一个像素）。如果我们想要更准确的计算，我们应该考虑这个案例并对其进行处理。

Cython代码

在本例中，我们将使用Cython。代码在下一张图片中显示，它是一大块代码。它的某些部分可以改进，但速度已经快得多。

代码的编写方式使得大部分代码完全在C API中运行。这种策略大大加快了执行速度，但也让我们可以利用Cython的并行模块。我们将在CPU中的多个核心之间拆分作业。

from libc.math cimport sin, cos, pi, ceil, floor, pow
from libc.stdlib cimport abort, malloc, free
import numpy as np
cimport numpy as np
cimport cython
from cython.parallel import prange, parallel
cimport openmp
 
 
cdef double get_pixel2d(
 double *image,
 Py_ssize_t n_rows, 
 Py_ssize_t n_cols,
 long x,
 long y) nogil:
 
 if (y < 0) or (y >= n_rows) or (x < 0) or (x >= n_cols):
 return 0
 else:
 return image[y * n_cols + x]
 
 
cdef double bilinear_interpolation(
 double *image,
 Py_ssize_t n_rows,
 Py_ssize_t n_cols,
 double x,
 double y) nogil:
 
 cdef double d_y, d_x, top_left, top_right, bottom_left, bottom_right
 cdef long min_y, min_x, max_y, max_x
 
 min_y = <long>floor(y)
 min_x = <long>floor(x)
 max_y = <long>ceil(y)
 max_x = <long>ceil(x)
 
 d_y = y - min_y
 d_x = x - min_x
 
 top_left = get_pixel2d(image, n_rows, n_cols, min_x, min_y)
 top_right = get_pixel2d(image, n_rows, n_cols, max_x, min_y)
 bottom_left = get_pixel2d(image, n_rows, n_cols, min_x, max_y)
 bottom_right = get_pixel2d(image, n_rows, n_cols, max_x, max_y)
 
 top = (1 - d_x) * top_left + d_x * top_right
 bottom = (1 - d_x) * bottom_left + d_x * bottom_right
 
 return (1 - d_y) * top + d_y * bottom
 
 
cdef double *joint_difference_distribution(
 double *image,
 Py_ssize_t n_rows,
 Py_ssize_t n_cols,
 int x0,
 int y0,
 int P,
 int R
) nogil:
 cdef Py_ssize_t p
 cdef double *T = <double *> malloc(sizeof(double) * P)
 cdef double x, y, gp, gc
 
 if T is NULL:
 abort()
 
 gc = get_pixel2d(image, n_rows, n_cols, x0, y0)
 
 for p in range(P):
 x = x0 + R * cos(2 * pi * p / P)
 y = y0 - R * sin(2 * pi * p / P)
 gp = bilinear_interpolation(image, n_rows, n_cols, x, y)
 T[p] = gp - gc
 
 return T
 
 
cdef int *binary_joint_distribution(double *T, Py_ssize_t T_size) nogil:
 cdef int *s_T = <int *> malloc(sizeof(int) * T_size)
 cdef Py_ssize_t i = 0
 
 for t in range(T_size):
 if T[t] >= 0.0:
 s_T[t] = 1
 else:
 s_T[t] = 0
 
 return s_T
 
 
cdef long LBP(double *T, int *s_T, Py_ssize_t T_size) nogil:
 cdef long LBP_pr = 0
 cdef Py_ssize_t i = 0
 
 for i in range(0, T_size):
 LBP_pr = LBP_pr + 2 ** i * s_T[i]
 
 return LBP_pr
 
 
cdef int is_uniform_pattern(int *s_T, Py_ssize_t s_T_size) nogil:
 cdef Py_ssize_t i = 0
 cdef int counter = 0
 
 for i in range(s_T_size - 1):
 if s_T[i] != s_T[i + 1]:
 counter += 1
 
 if counter > 2:
 return 0
 return 1
 
 
cdef int create_index(int *s_T, Py_ssize_t s_T_size) nogil:
 cdef int n_ones = 0
 cdef int rot_index = -1
 cdef int first_one = -1
 cdef int first_zero = -1
 cdef int lbp = -1
 
 cdef Py_ssize_t i
 for i in range(s_T_size):
 if s_T[i]:
 n_ones += 1
 if first_one == -1:
 first_one = i
 else:
 if first_zero == -1:
 first_zero = i
 
 if n_ones == 0:
 lbp = 0
 elif n_ones == s_T_size:
 lbp = s_T_size * (s_T_size - 1) + 1
 else:
 if first_one == 0:
 rot_index = n_ones - first_zero
 else:
 rot_index = s_T_size - first_one
 lbp = 1 + (n_ones - 1) * s_T_size + rot_index
 return lbp
 
 
cdef int LBP_uniform(int *s_T, Py_ssize_t s_T_size) nogil:
 cdef int LBP_pru = 0
 cdef Py_ssize_t i = 0
 
 if is_uniform_pattern(s_T, s_T_size):
 LBP_pru = create_index(s_T, s_T_size)
 else:
 LBP_pru = 2 + s_T_size * (s_T_size - 1)
 
 return LBP_pru
 
 
@cython.boundscheck(False)
@cython.wraparound(False)
def local_binary_patterns(
 double[:, ::1] image,
 int P,
 int R,
 int num_threads=1
):
 
 cdef Py_ssize_t x = 0
 cdef Py_ssize_t y = 0
 cdef int n_rows = image.shape[0]
 cdef int n_cols = image.shape[1]
 cdef int[:, ::1] lbp = np.zeros([n_rows, n_cols], dtype=np.int32) 
 
 with nogil, parallel(num_threads=num_threads):
 for y in prange(n_rows, schedule='static'):
 for x in prange(n_cols, schedule='static'):
 T = joint_difference_distribution(&image[0][0], n_rows, n_cols, x, y, P, R)
 s_T = binary_joint_distribution(T, P)
 lbp[y, x] = LBP_uniform(s_T, P)
 
 return np.asarray(lbp)

使用4个线程，我们可以在不到150毫秒的时间内计算所有像素的局部二进制模式。

与相似图像的比较

让我们拍摄另一张砖块图片，但这张图片会有不同的纹理。

两个直方图非常相似，它们应该是，最后它们都是砖块。尽管如此，20到40的特征在两个图像中都非常不同。这意味着通过良好的机器学习算法，我们可以正确地对它们进行分类。

结论

局部二进制模式是简单但有效的机器学习特征。背后的理论并不难理解，并且易于编码。然而，如果我们完全使用Python编写代码，我们将遇到一些性能问题。我们用Cython解决了这个问题，结果非常令人印象深刻。