局部二值模式(LBP)的matlab实现

LBP是一种典型的图像局部纹理特征描述子，其相关原理网上已经有不少不错的博客介绍过了，如https://blog.csdn.net/xpplearnc/article/details/40658211，这里本人就不再赘述了。下面给出本人的matlab实现及示例结果以作示范（第二幅测试图为随手截取的网上的图片，由于半径及邻域点数不同、是否对均匀模式归化及归化到何值方面可能存在分歧，所以结果并不完全一致，这里只是做个示范。附图也不做介绍了，大家自己明白就好）：

function flag=isUniformMode(x,n)
flag=1; %是否均匀模式
S=zeros(n+1,1);
m=x;  %余数
for j=n-1:-1:0 %对比特位的循环放在外层可以提升效率，这里这样写是为可读性起见
    if m<2^j
        S(j+2)=0;
    else
        S(j+2)=1;
        m=m-2^j;
    end
end
S(1)=S(n+1);
num=0;
for j=1:n
    if S(j+1)~=S(j)
        num=num+1;
    end
    if num>2
        flag=0;
        break;
    end
end

function map=mapOfP2shiftedP(n)
%构建由移位前LBP值到最小移位LBP值的map
Key=num2cell(0:2^n-1); 
matVal=zeros(2^n,1);
NoneUniform=0;
for i=0:n-1
    isEven=mod(i,2);
    NoneUniform=NoneUniform+isEven*2^i; %构造一个非均匀模式
end
for i=0:2^n-1
    val=zeros(n,1);
    val(n)=i;
    % 求移位后的取值，断点从高位往低位移动
    for j=n-1:-1:1
        p=0; %高位部分
        q=0; %低位部分
        if i<2^j
            q=i;
            p=0;
        else
            q=mod(i,2^j);
            p=(i-q)/2^j;
        end
        val(j)=q*2^(n-j)+p;
    end
    minVal=min(val);
    if isUniformMode(minVal,n) 
        matVal(i+1)=minVal;
    else
        matVal(i+1)=NoneUniform; %非均匀模式统一为一种模式
    end
end
Val=num2cell(matVal);
map=containers.Map(Key,Val);

function out=LBP(img,radius,n)
%构建LBP图像
%使用圆形邻域，等分选取n点

%构造移位前后lbp映射表
map=mapOfP2shiftedP(n);

%构建相对位置与权重阵列
dx1=zeros(n,1);
dy1=zeros(n,1);
dx2=zeros(n,1);
dy2=zeros(n,1);
w11=zeros(n,1);
w12=zeros(n,1);
w21=zeros(n,1);
w22=zeros(n,1);
for i=0:n-1
    x=radius*cos(i*2*pi/n);
    y=radius*sin(i*2*pi/n);
    dx1(i+1)=floor(x);
    dx2(i+1)=ceil(x);
    dy1(i+1)=floor(y);
    dy2(i+1)=ceil(y);
    wx=ceil(x)-x;
    wy=ceil(y)-y;
    w11(i+1)=wx*wy;
    w21(i+1)=(1-wx)*wy;
    w12(i+1)=wx*(1-wy);
    w22(i+1)=(1-wx)*(1-wy);
end

out=zeros(size(img,1)-2*radius,size(img,2)-2*radius);
for y=radius+1:size(img,1)-radius
    for x=radius+1:size(img,2)-radius
        curVal=img(y,x);
        num=0;
        for i=1:n
            intVal=w11(i)*img(y+dy1(i),x+dx1(i))+w21(i)*img(y+dy1(i),x+dx2(i))+...
                       w12(i)*img(y+dy2(i),x+dx1(i))+w22(i)*img(y+dy2(i),x+dx2(i));
            if intVal>curVal
                num=num+2^(i-1);
            end
        end
        out(y-radius,x-radius)=map(num);
    end
end

function h=LBPhist(LBP,n)
%获得LBP统计直方图

%构造移位前后lbp映射表
map=mapOfP2shiftedP(n);

%构造移位后lbp到直方图bin的映射
v=cell2mat(values(map));
v=unique(v);
key=num2cell(v);
value=num2cell(1:length(v));
map2=containers.Map(key,value);

h=zeros(length(v),1);
%直方图统计
for i=1:size(LBP,1)
    for j=1:size(LBP,2)
        idx=map2(LBP(i,j));
        h(idx)=h(idx)+1;
    end
end

 脚本文件：

clear all;close all;clc;
img=rgb2gray(imread('lena.tiff'));
radius=2;
n=16;
lbp=LBP(img,radius,n);
figure,imshow(lbp/(2^n-1));
imwrite(lbp/(2^n-1),'res.bmp');
h=LBPhist(lbp,n);
figure,bar(h);
title('LBP直方图');

 测试图及运行结果：