Local Binary Pattern for texture classificationΒΆ

In this example, we will see how to classify textures based on LBP (Local Binary Pattern). LBP looks at points surrounding a central point and tests whether the surrounding points are greater than or less than the central point (i.e. gives a binary result).

Before trying out LBP on an image, it helps to look at a schematic of LBPs. The below code is just used to plot the schematic.

from __future__ import print_function
import numpy as np
import matplotlib.pyplot as plt


METHOD = 'uniform'
plt.rcParams['font.size'] = 9


def plot_circle(ax, center, radius, color):
    circle = plt.Circle(center, radius, facecolor=color, edgecolor='0.5')
    ax.add_patch(circle)


def plot_lbp_model(ax, binary_values):
    """Draw the schematic for a local binary pattern."""
    # Geometry spec
    theta = np.deg2rad(45)
    R = 1
    r = 0.15
    w = 1.5
    gray = '0.5'

    # Draw the central pixel.
    plot_circle(ax, (0, 0), radius=r, color=gray)
    # Draw the surrounding pixels.
    for i, facecolor in enumerate(binary_values):
        x = R * np.cos(i * theta)
        y = R * np.sin(i * theta)
        plot_circle(ax, (x, y), radius=r, color=str(facecolor))

    # Draw the pixel grid.
    for x in np.linspace(-w, w, 4):
        ax.axvline(x, color=gray)
        ax.axhline(x, color=gray)

    # Tweak the layout.
    ax.axis('image')
    ax.axis('off')
    size = w + 0.2
    ax.set_xlim(-size, size)
    ax.set_ylim(-size, size)


fig, axes = plt.subplots(ncols=5, figsize=(7, 2))

titles = ['flat', 'flat', 'edge', 'corner', 'non-uniform']

binary_patterns = [np.zeros(8),
                   np.ones(8),
                   np.hstack([np.ones(4), np.zeros(4)]),
                   np.hstack([np.zeros(3), np.ones(5)]),
                   [1, 0, 0, 1, 1, 1, 0, 0]]

for ax, values, name in zip(axes, binary_patterns, titles):
    plot_lbp_model(ax, values)
    ax.set_title(name)

../_images/plot_local_binary_pattern_1.png

The figure above shows example results with black (or white) representing pixels that are less (or more) intense than the central pixel. When surrounding pixels are all black or all white, then that image region is flat (i.e. featureless). Groups of continuous black or white pixels are considered “uniform” patterns that can be interpreted as corners or edges. If pixels switch back-and-forth between black and white pixels, the pattern is considered “non-uniform”.

When using LBP to detect texture, you measure a collection of LBPs over an image patch and look at the distribution of these LBPs. Lets apply LBP to a brick texture.

from skimage.transform import rotate
from skimage.feature import local_binary_pattern
from skimage import data
from skimage.color import label2rgb

# settings for LBP
radius = 3
n_points = 8 * radius


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)


def highlight_bars(bars, indexes):
    for i in indexes:
        bars[i].set_facecolor('r')


image = data.load('brick.png')
lbp = local_binary_pattern(image, n_points, radius, METHOD)

def hist(ax, lbp):
    n_bins = lbp.max() + 1
    return ax.hist(lbp.ravel(), normed=True, bins=n_bins, range=(0, n_bins),
                   facecolor='0.5')

# plot histograms of LBP of textures
fig, (ax_img, ax_hist) = plt.subplots(nrows=2, ncols=3, figsize=(9, 6))
plt.gray()

titles = ('edge', 'flat', 'corner')
w = width = radius - 1
edge_labels = range(n_points // 2 - w, n_points // 2 + w + 1)
flat_labels = list(range(0, w + 1)) + list(range(n_points - w, n_points + 2))
i_14 = n_points // 4            # 1/4th of the histogram
i_34 = 3 * (n_points // 4)      # 3/4th of the histogram
corner_labels = (list(range(i_14 - w, i_14 + w + 1)) +
                 list(range(i_34 - w, i_34 + w + 1)))

label_sets = (edge_labels, flat_labels, corner_labels)

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(ymax=np.max(counts[:-1]))
    ax.set_xlim(xmax=n_points + 2)
    ax.set_title(name)

ax_hist[0].set_ylabel('Percentage')
for ax in ax_img:
    ax.axis('off')

../_images/plot_local_binary_pattern_2.png

The above plot highlights flat, edge-like, and corner-like regions of the image.

The histogram of the LBP result is a good measure to classify textures. Here, we test the histogram distributions against each other using the Kullback-Leibler-Divergence.

# settings for 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 = lbp.max() + 1
    hist, _ = np.histogram(lbp, normed=True, bins=n_bins, range=(0, n_bins))
    for name, ref in refs.items():
        ref_hist, _ = np.histogram(ref, normed=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.load('brick.png')
grass = data.load('grass.png')
wall = data.load('rough-wall.png')

refs = {
    'brick': local_binary_pattern(brick, n_points, radius, METHOD),
    'grass': local_binary_pattern(grass, n_points, radius, METHOD),
    'wall': local_binary_pattern(wall, 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(wall)
ax3.axis('off')
hist(ax6, refs['wall'])

../_images/plot_local_binary_pattern_3.png

plt.show()

STDOUT


        

STDERR


        

Python source code: download (generated using skimage 0.11dev)

IPython Notebook: download (generated using skimage 0.11dev)

ZnJvbSBfX2Z1dHVyZV9fIGltcG9ydCBwcmludF9mdW5jdGlvbgppbXBvcnQgbnVtcHkgYXMgbnAKaW1wb3J0IG1hdHBsb3RsaWIucHlwbG90IGFzIHBsdAoKCk1FVEhPRCA9ICd1bmlmb3JtJwpwbHQucmNQYXJhbXNbJ2ZvbnQuc2l6ZSddID0gOQoKCmRlZiBwbG90X2NpcmNsZShheCwgY2VudGVyLCByYWRpdXMsIGNvbG9yKToKICAgIGNpcmNsZSA9IHBsdC5DaXJjbGUoY2VudGVyLCByYWRpdXMsIGZhY2Vjb2xvcj1jb2xvciwgZWRnZWNvbG9yPScwLjUnKQogICAgYXguYWRkX3BhdGNoKGNpcmNsZSkKCgpkZWYgcGxvdF9sYnBfbW9kZWwoYXgsIGJpbmFyeV92YWx1ZXMpOgogICAgIiIiRHJhdyB0aGUgc2NoZW1hdGljIGZvciBhIGxvY2FsIGJpbmFyeSBwYXR0ZXJuLiIiIgogICAgIyBHZW9tZXRyeSBzcGVjCiAgICB0aGV0YSA9IG5wLmRlZzJyYWQoNDUpCiAgICBSID0gMQogICAgciA9IDAuMTUKICAgIHcgPSAxLjUKICAgIGdyYXkgPSAnMC41JwoKICAgICMgRHJhdyB0aGUgY2VudHJhbCBwaXhlbC4KICAgIHBsb3RfY2lyY2xlKGF4LCAoMCwgMCksIHJhZGl1cz1yLCBjb2xvcj1ncmF5KQogICAgIyBEcmF3IHRoZSBzdXJyb3VuZGluZyBwaXhlbHMuCiAgICBmb3IgaSwgZmFjZWNvbG9yIGluIGVudW1lcmF0ZShiaW5hcnlfdmFsdWVzKToKICAgICAgICB4ID0gUiAqIG5wLmNvcyhpICogdGhldGEpCiAgICAgICAgeSA9IFIgKiBucC5zaW4oaSAqIHRoZXRhKQogICAgICAgIHBsb3RfY2lyY2xlKGF4LCAoeCwgeSksIHJhZGl1cz1yLCBjb2xvcj1zdHIoZmFjZWNvbG9yKSkKCiAgICAjIERyYXcgdGhlIHBpeGVsIGdyaWQuCiAgICBmb3IgeCBpbiBucC5saW5zcGFjZSgtdywgdywgNCk6CiAgICAgICAgYXguYXh2bGluZSh4LCBjb2xvcj1ncmF5KQogICAgICAgIGF4LmF4aGxpbmUoeCwgY29sb3I9Z3JheSkKCiAgICAjIFR3ZWFrIHRoZSBsYXlvdXQuCiAgICBheC5heGlzKCdpbWFnZScpCiAgICBheC5heGlzKCdvZmYnKQogICAgc2l6ZSA9IHcgKyAwLjIKICAgIGF4LnNldF94bGltKC1zaXplLCBzaXplKQogICAgYXguc2V0X3lsaW0oLXNpemUsIHNpemUpCgoKZmlnLCBheGVzID0gcGx0LnN1YnBsb3RzKG5jb2xzPTUsIGZpZ3NpemU9KDcsIDIpKQoKdGl0bGVzID0gWydmbGF0JywgJ2ZsYXQnLCAnZWRnZScsICdjb3JuZXInLCAnbm9uLXVuaWZvcm0nXQoKYmluYXJ5X3BhdHRlcm5zID0gW25wLnplcm9zKDgpLAogICAgICAgICAgICAgICAgICAgbnAub25lcyg4KSwKICAgICAgICAgICAgICAgICAgIG5wLmhzdGFjayhbbnAub25lcyg0KSwgbnAuemVyb3MoNCldKSwKICAgICAgICAgICAgICAgICAgIG5wLmhzdGFjayhbbnAuemVyb3MoMyksIG5wLm9uZXMoNSldKSwKICAgICAgICAgICAgICAgICAgIFsxLCAwLCAwLCAxLCAxLCAxLCAwLCAwXV0KCmZvciBheCwgdmFsdWVzLCBuYW1lIGluIHppcChheGVzLCBiaW5hcnlfcGF0dGVybnMsIHRpdGxlcyk6CiAgICBwbG90X2xicF9tb2RlbChheCwgdmFsdWVzKQogICAgYXguc2V0X3RpdGxlKG5hbWUp
CmZyb20gc2tpbWFnZS50cmFuc2Zvcm0gaW1wb3J0IHJvdGF0ZQpmcm9tIHNraW1hZ2UuZmVhdHVyZSBpbXBvcnQgbG9jYWxfYmluYXJ5X3BhdHRlcm4KZnJvbSBza2ltYWdlIGltcG9ydCBkYXRhCmZyb20gc2tpbWFnZS5jb2xvciBpbXBvcnQgbGFiZWwycmdiCgojIHNldHRpbmdzIGZvciBMQlAKcmFkaXVzID0gMwpuX3BvaW50cyA9IDggKiByYWRpdXMKCgpkZWYgb3ZlcmxheV9sYWJlbHMoaW1hZ2UsIGxicCwgbGFiZWxzKToKICAgIG1hc2sgPSBucC5sb2dpY2FsX29yLnJlZHVjZShbbGJwID09IGVhY2ggZm9yIGVhY2ggaW4gbGFiZWxzXSkKICAgIHJldHVybiBsYWJlbDJyZ2IobWFzaywgaW1hZ2U9aW1hZ2UsIGJnX2xhYmVsPTAsIGFscGhhPTAuNSkKCgpkZWYgaGlnaGxpZ2h0X2JhcnMoYmFycywgaW5kZXhlcyk6CiAgICBmb3IgaSBpbiBpbmRleGVzOgogICAgICAgIGJhcnNbaV0uc2V0X2ZhY2Vjb2xvcigncicpCgoKaW1hZ2UgPSBkYXRhLmxvYWQoJ2JyaWNrLnBuZycpCmxicCA9IGxvY2FsX2JpbmFyeV9wYXR0ZXJuKGltYWdlLCBuX3BvaW50cywgcmFkaXVzLCBNRVRIT0QpCgpkZWYgaGlzdChheCwgbGJwKToKICAgIG5fYmlucyA9IGxicC5tYXgoKSArIDEKICAgIHJldHVybiBheC5oaXN0KGxicC5yYXZlbCgpLCBub3JtZWQ9VHJ1ZSwgYmlucz1uX2JpbnMsIHJhbmdlPSgwLCBuX2JpbnMpLAogICAgICAgICAgICAgICAgICAgZmFjZWNvbG9yPScwLjUnKQoKIyBwbG90IGhpc3RvZ3JhbXMgb2YgTEJQIG9mIHRleHR1cmVzCmZpZywgKGF4X2ltZywgYXhfaGlzdCkgPSBwbHQuc3VicGxvdHMobnJvd3M9MiwgbmNvbHM9MywgZmlnc2l6ZT0oOSwgNikpCnBsdC5ncmF5KCkKCnRpdGxlcyA9ICgnZWRnZScsICdmbGF0JywgJ2Nvcm5lcicpCncgPSB3aWR0aCA9IHJhZGl1cyAtIDEKZWRnZV9sYWJlbHMgPSByYW5nZShuX3BvaW50cyAvLyAyIC0gdywgbl9wb2ludHMgLy8gMiArIHcgKyAxKQpmbGF0X2xhYmVscyA9IGxpc3QocmFuZ2UoMCwgdyArIDEpKSArIGxpc3QocmFuZ2Uobl9wb2ludHMgLSB3LCBuX3BvaW50cyArIDIpKQppXzE0ID0gbl9wb2ludHMgLy8gNCAgICAgICAgICAgICMgMS80dGggb2YgdGhlIGhpc3RvZ3JhbQppXzM0ID0gMyAqIChuX3BvaW50cyAvLyA0KSAgICAgICMgMy80dGggb2YgdGhlIGhpc3RvZ3JhbQpjb3JuZXJfbGFiZWxzID0gKGxpc3QocmFuZ2UoaV8xNCAtIHcsIGlfMTQgKyB3ICsgMSkpICsKICAgICAgICAgICAgICAgICBsaXN0KHJhbmdlKGlfMzQgLSB3LCBpXzM0ICsgdyArIDEpKSkKCmxhYmVsX3NldHMgPSAoZWRnZV9sYWJlbHMsIGZsYXRfbGFiZWxzLCBjb3JuZXJfbGFiZWxzKQoKZm9yIGF4LCBsYWJlbHMgaW4gemlwKGF4X2ltZywgbGFiZWxfc2V0cyk6CiAgICBheC5pbXNob3cob3ZlcmxheV9sYWJlbHMoaW1hZ2UsIGxicCwgbGFiZWxzKSkKCmZvciBheCwgbGFiZWxzLCBuYW1lIGluIHppcChheF9oaXN0LCBsYWJlbF9zZXRzLCB0aXRsZXMpOgogICAgY291bnRzLCBfLCBiYXJzID0gaGlzdChheCwgbGJwKQogICAgaGlnaGxpZ2h0X2JhcnMoYmFycywgbGFiZWxzKQogICAgYXguc2V0X3lsaW0oeW1heD1ucC5tYXgoY291bnRzWzotMV0pKQogICAgYXguc2V0X3hsaW0oeG1heD1uX3BvaW50cyArIDIpCiAgICBheC5zZXRfdGl0bGUobmFtZSkKCmF4X2hpc3RbMF0uc2V0X3lsYWJlbCgnUGVyY2VudGFnZScpCmZvciBheCBpbiBheF9pbWc6CiAgICBheC5heGlzKCdvZmYnKQ==
CiMgc2V0dGluZ3MgZm9yIExCUApyYWRpdXMgPSAyCm5fcG9pbnRzID0gOCAqIHJhZGl1cwoKCmRlZiBrdWxsYmFja19sZWlibGVyX2RpdmVyZ2VuY2UocCwgcSk6CiAgICBwID0gbnAuYXNhcnJheShwKQogICAgcSA9IG5wLmFzYXJyYXkocSkKICAgIGZpbHQgPSBucC5sb2dpY2FsX2FuZChwICE9IDAsIHEgIT0gMCkKICAgIHJldHVybiBucC5zdW0ocFtmaWx0XSAqIG5wLmxvZzIocFtmaWx0XSAvIHFbZmlsdF0pKQoKCmRlZiBtYXRjaChyZWZzLCBpbWcpOgogICAgYmVzdF9zY29yZSA9IDEwCiAgICBiZXN0X25hbWUgPSBOb25lCiAgICBsYnAgPSBsb2NhbF9iaW5hcnlfcGF0dGVybihpbWcsIG5fcG9pbnRzLCByYWRpdXMsIE1FVEhPRCkKICAgIG5fYmlucyA9IGxicC5tYXgoKSArIDEKICAgIGhpc3QsIF8gPSBucC5oaXN0b2dyYW0obGJwLCBub3JtZWQ9VHJ1ZSwgYmlucz1uX2JpbnMsIHJhbmdlPSgwLCBuX2JpbnMpKQogICAgZm9yIG5hbWUsIHJlZiBpbiByZWZzLml0ZW1zKCk6CiAgICAgICAgcmVmX2hpc3QsIF8gPSBucC5oaXN0b2dyYW0ocmVmLCBub3JtZWQ9VHJ1ZSwgYmlucz1uX2JpbnMsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgcmFuZ2U9KDAsIG5fYmlucykpCiAgICAgICAgc2NvcmUgPSBrdWxsYmFja19sZWlibGVyX2RpdmVyZ2VuY2UoaGlzdCwgcmVmX2hpc3QpCiAgICAgICAgaWYgc2NvcmUgPCBiZXN0X3Njb3JlOgogICAgICAgICAgICBiZXN0X3Njb3JlID0gc2NvcmUKICAgICAgICAgICAgYmVzdF9uYW1lID0gbmFtZQogICAgcmV0dXJuIGJlc3RfbmFtZQoKCmJyaWNrID0gZGF0YS5sb2FkKCdicmljay5wbmcnKQpncmFzcyA9IGRhdGEubG9hZCgnZ3Jhc3MucG5nJykKd2FsbCA9IGRhdGEubG9hZCgncm91Z2gtd2FsbC5wbmcnKQoKcmVmcyA9IHsKICAgICdicmljayc6IGxvY2FsX2JpbmFyeV9wYXR0ZXJuKGJyaWNrLCBuX3BvaW50cywgcmFkaXVzLCBNRVRIT0QpLAogICAgJ2dyYXNzJzogbG9jYWxfYmluYXJ5X3BhdHRlcm4oZ3Jhc3MsIG5fcG9pbnRzLCByYWRpdXMsIE1FVEhPRCksCiAgICAnd2FsbCc6IGxvY2FsX2JpbmFyeV9wYXR0ZXJuKHdhbGwsIG5fcG9pbnRzLCByYWRpdXMsIE1FVEhPRCkKfQoKIyBjbGFzc2lmeSByb3RhdGVkIHRleHR1cmVzCnByaW50KCdSb3RhdGVkIGltYWdlcyBtYXRjaGVkIGFnYWluc3QgcmVmZXJlbmNlcyB1c2luZyBMQlA6JykKcHJpbnQoJ29yaWdpbmFsOiBicmljaywgcm90YXRlZDogMzBkZWcsIG1hdGNoIHJlc3VsdDogJywKICAgICAgbWF0Y2gocmVmcywgcm90YXRlKGJyaWNrLCBhbmdsZT0zMCwgcmVzaXplPUZhbHNlKSkpCnByaW50KCdvcmlnaW5hbDogYnJpY2ssIHJvdGF0ZWQ6IDcwZGVnLCBtYXRjaCByZXN1bHQ6ICcsCiAgICAgIG1hdGNoKHJlZnMsIHJvdGF0ZShicmljaywgYW5nbGU9NzAsIHJlc2l6ZT1GYWxzZSkpKQpwcmludCgnb3JpZ2luYWw6IGdyYXNzLCByb3RhdGVkOiAxNDVkZWcsIG1hdGNoIHJlc3VsdDogJywKICAgICAgbWF0Y2gocmVmcywgcm90YXRlKGdyYXNzLCBhbmdsZT0xNDUsIHJlc2l6ZT1GYWxzZSkpKQoKIyBwbG90IGhpc3RvZ3JhbXMgb2YgTEJQIG9mIHRleHR1cmVzCmZpZywgKChheDEsIGF4MiwgYXgzKSwgKGF4NCwgYXg1LCBheDYpKSA9IHBsdC5zdWJwbG90cyhucm93cz0yLCBuY29scz0zLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgZmlnc2l6ZT0oOSwgNikpCnBsdC5ncmF5KCkKCmF4MS5pbXNob3coYnJpY2spCmF4MS5heGlzKCdvZmYnKQpoaXN0KGF4NCwgcmVmc1snYnJpY2snXSkKYXg0LnNldF95bGFiZWwoJ1BlcmNlbnRhZ2UnKQoKYXgyLmltc2hvdyhncmFzcykKYXgyLmF4aXMoJ29mZicpCmhpc3QoYXg1LCByZWZzWydncmFzcyddKQpheDUuc2V0X3hsYWJlbCgnVW5pZm9ybSBMQlAgdmFsdWVzJykKCmF4My5pbXNob3cod2FsbCkKYXgzLmF4aXMoJ29mZicpCmhpc3QoYXg2LCByZWZzWyd3YWxsJ10p
CnBsdC5zaG93KCk=