LBP源代码阅读

静下心来看了一下LBP的源代码，作者写得很漂亮。以前经常用，但没有注意过细节，写了一些注释与大家分享。

%LBP returns the local binary pattern image or LBP histogram of an image.%  J = LBP(I,R,N,MAPPING,MODE) returns either a local binary pattern%  coded image or the local binary pattern histogram of an intensity%  image I. The LBP codes are computed using N sampling points on a %  circle of radius R and using mapping table defined by MAPPING. %  See the getmapping function for different mappings and use 0 for%  no mapping. Possible values for MODE are%       'h' or 'hist'  to get a histogram of LBP codes%       'nh'           to get a normalized histogram%  Otherwise an LBP code image is returned.%%  J = LBP(I) returns the original (basic) LBP histogram of image I%%  J = LBP(I,SP,MAPPING,MODE) computes the LBP codes using n sampling%  points defined in (n * 2) matrix SP. The sampling points should be%  defined around the origin (coordinates (0,0)).%%  Examples%  --------%       I=imread('test1.bmp');%       mapping=getmapping(8,'u2'); %       H1=lbp(I,1,8,mapping,'h'); %LBP histogram in (8,1) neighborhood%                                  %using uniform patterns%       subplot(2,1,1),stem(H1);% %       H2=lbp(I);%       subplot(2,1,2),stem(H2);% %       SP=[-1 -1; -1 0; -1 1; 0 -1; -0 1; 1 -1; 1 0; 1 1];%       I2=lbp(I,SP,0,'i'); %LBP code image using sampling points in SP%                           %and no mapping. Now H2 is equal to histogram%                           %of I2.function result = lbp(varargin) % image,radius,neighbors,mapping,mode)% Version 0.3.2% Authors: Marko Heikkil?and Timo Ahonen% Changelog% Version 0.3.2: A bug fix to enable using mappings together with a% predefined spoints array% Version 0.3.1: Changed MAPPING input to be a struct containing the mapping% table and the number of bins to make the function run faster with high number% of sampling points. Lauge Sorensen is acknowledged for spotting this problem.% Check number of input arguments.% %% 检查参数的个数nargin，使其大于1小于5。如果不在此区间，就报错error(nargchk(1,5,nargin));image=varargin{1};% % %% 把第一个参数赋值给imaged_image=double(image);% % % 把图像从uint8转成double类型，以便以后计算% % % 只有给出待处理的图像（一个参数）时，使用默认的设置。% % % sp定义了中心点与它的近邻的相对位置% % % neighbors定义近邻个数% % % mapping定义的映射% % % mode区别直方图的类型，'h' or 'hist'是直方图，nh是规一化的直方图if nargin==1    spoints=[-1 -1; -1 0; -1 1; 0 -1; -0 1; 1 -1; 1 0; 1 1];                   neighbors=8;    mapping=0;    mode='h';end% % % 给出两个参数，并且第二个参数（代表近邻半径）的长度为1时% % % 只给出了近邻的半径，没给出近邻的个数，报错if (nargin == 2) && (length(varargin{2}) == 1)    error('Input arguments');end% % % 如果给出两个以上的参数，并且第二个参数（代表近邻半径）的长度为1% % % 半径设为第二个参数% % % 近邻个数设为第三个参数if (nargin > 2) && (length(varargin{2}) == 1)    radius=varargin{2};    neighbors=varargin{3};        spoints=zeros(neighbors,2);    % % %     把360度均匀分成neighbors分，以计算近邻点与中心的相对坐标    % Angle step.    a = 2*pi/neighbors;    % % %       计算坐标，每一维代表y,第二维代表x% % %       spoints的第i行代表第i个近邻    for i = 1:neighbors        spoints(i,1) = -radius*sin((i-1)*a);        spoints(i,2) = radius*cos((i-1)*a);    end% % %     如果参数个数大于等于4，第四个参数赋值给映射mapping；否则，无映射。    if(nargin >= 4)        mapping=varargin{4};        if(isstruct(mapping) && mapping.samples ~= neighbors)            error('Incompatible mapping');        end    else        mapping=0;    end% % %     第五个参数确定直方图的属性    if(nargin >= 5)        mode=varargin{5};    else        mode='h';    endend% % % 如果参数个数大于1，并且第二个参数的长度大于1。则第二个参数给出近邻点与中心点的相对位置if (nargin > 1) && (length(varargin{2}) > 1)    spoints=varargin{2};    neighbors=size(spoints,1);% % %     如果还有第三个参数，把它赋值给映射mapping    if(nargin >= 3)        mapping=varargin{3};        if(isstruct(mapping) && mapping.samples ~= neighbors)            error('Incompatible mapping');        end    else        mapping=0;    end        if(nargin >= 4)        mode=varargin{4};    else        mode='h';    end   end% Determine the dimensions of the input image.图像的大小，第一维是y，第二维是x[ysize xsize] = size(image);% % % 确定block的左上和右下两个点miny=min(spoints(:,1));maxy=max(spoints(:,1));minx=min(spoints(:,2));maxx=max(spoints(:,2));% Block size, each LBP code is computed within a block of size% bsizey*bsizex % % % block的大小bsizey=ceil(max(maxy,0))-floor(min(miny,0))+1;bsizex=ceil(max(maxx,0))-floor(min(minx,0))+1;% Coordinates of origin (0,0) in the block% % % 在block里中心点的坐标origy=1-floor(min(miny,0));origx=1-floor(min(minx,0));% Minimum allowed size for the input image depends% on the radius of the used LBP operator.% % % 检查block和img的大小if(xsize < bsizex || ysize < bsizey)  error('Too small input image. Should be at least (2*radius+1) x (2*radius+1)');end% Calculate dx and dy;dx = xsize - bsizex;dy = ysize - bsizey;% Fill the center pixel matrix C.% % % 所有可以作为模板中心点的像素集合C = image(origy:origy+dy,origx:origx+dx);d_C = double(C);bins = 2^neighbors;% Initialize the result matrix with zeros.result=zeros(dy+1,dx+1);% % % 初始化结果矩阵%Compute the LBP code image  这一段写得很漂亮！！！！% % % 对于每一个neighbor，先使要比较的点与中心点对齐，然后利用D = N >= C比较它们的大小。for i = 1:neighbors  y = spoints(i,1)+origy;  x = spoints(i,2)+origx;  % Calculate floors, ceils and rounds for the x and y.  fy = floor(y); cy = ceil(y); ry = round(y);  fx = floor(x); cx = ceil(x); rx = round(x);  % Check if interpolation is needed.  if (abs(x - rx) < 1e-6) && (abs(y - ry) < 1e-6)    % Interpolation is not needed, use original datatypes    N = image(ry:ry+dy,rx:rx+dx);    D = N >= C;   else    % Interpolation needed, use double type images     ty = y - fy;    tx = x - fx;    % Calculate the interpolation weights.    w1 = (1 - tx) * (1 - ty);    w2 =      tx  * (1 - ty);    w3 = (1 - tx) *      ty ;    w4 =      tx  *      ty ;    % Compute interpolated pixel values    N = w1*d_image(fy:fy+dy,fx:fx+dx) + w2*d_image(fy:fy+dy,cx:cx+dx) + ...        w3*d_image(cy:cy+dy,fx:fx+dx) + w4*d_image(cy:cy+dy,cx:cx+dx);    D = N >= d_C;   end    % Update the result matrix.% % %   更新结果矩阵  v = 2^(i-1);  result = result + v*D;end%Apply mapping if it is defined% % % 如果mapping已经存在，那么利用这个mapping.if isstruct(mapping)    bins = mapping.num;    for i = 1:size(result,1)        for j = 1:size(result,2)            result(i,j) = mapping.table(result(i,j)+1);        end    endend% % % 如果要参数列表指定了直方图的属性，计算直方图if (strcmp(mode,'h') || strcmp(mode,'hist') || strcmp(mode,'nh'))    % Return with LBP histogram if mode equals 'hist'.    result=hist(result(:),0:(bins-1));    if (strcmp(mode,'nh'))        result=result/sum(result);    endelse    %Otherwise return a matrix of unsigned integers% % %     如果没有指定直方图的属性，返回数值方阵    if ((bins-1)<=intmax('uint8'))        result=uint8(result);    elseif ((bins-1)<=intmax('uint16'))        result=uint16(result);    else        result=uint32(result);    endendend

生成mapping的函数

%GETMAPPING returns a structure containing a mapping table for LBP codes.%  MAPPING = GETMAPPING(SAMPLES,MAPPINGTYPE) returns a %  structure containing a mapping table for%  LBP codes in a neighbourhood of SAMPLES sampling%  points. Possible values for MAPPINGTYPE are%       'u2'   for uniform LBP%       'ri'   for rotation-invariant LBP%       'riu2' for uniform rotation-invariant LBP.%%  Example:%       I=imread('rice.tif');%       MAPPING=getmapping(16,'riu2');%       LBPHIST=lbp(I,2,16,MAPPING,'hist');%  Now LBPHIST contains a rotation-invariant uniform LBP%  histogram in a (16,2) neighbourhood.%function mapping = getmapping(samples,mappingtype)% Version 0.1.1% Authors: Marko Heikkil?and Timo Ahonen% Changelog% 0.1.1 Changed output to be a structure% Fixed a bug causing out of memory errors when generating rotation % invariant mappings with high number of sampling points.% Lauge Sorensen is acknowledged for spotting this problem.table = 0:2^samples-1;newMax  = 0; %number of patterns in the resulting LBP codeindex   = 0;if strcmp(mappingtype,'u2') %Uniform 2  newMax = samples*(samples-1) + 3;   for i = 0:2^samples-1% % %       j是i循环左移的结果    j = bitset(bitshift(i,1,'uint8'),1,bitget(i,samples)); %rotate left% % %     计算跳变的次数    numt = sum(bitget(bitxor(i,j),1:samples)); %number of 1->0 and                                               %0->1 transitions                                               %in binary string                                                %x is equal to the                                               %number of 1-bits in                                               %XOR(x,Rotate left(x)) % % %     如果跳变次数不大于2，那么新建一个标记index;否则放入最后一类    if numt <= 2      table(i+1) = index;      index = index + 1;    else      table(i+1) = newMax - 1;    end  endendif strcmp(mappingtype,'ri') %Rotation invariant  tmpMap = zeros(2^samples,1) - 1;  for i = 0:2^samples-1    rm = i;    r  = i;% % %     计算所有rotate中最小的一个    for j = 1:samples-1      r = bitset(bitshift(r,1,'uint8'),1,bitget(r,samples)); %rotate                                                             %left      if r < rm        rm = r;      end    end% % %     同上    if tmpMap(rm+1) < 0      tmpMap(rm+1) = newMax;      newMax = newMax + 1;    end    table(i+1) = tmpMap(rm+1);  endendif strcmp(mappingtype,'riu2') %Uniform & Rotation invariant  newMax = samples + 2;  for i = 0:2^samples - 1    j = bitset(bitshift(i,1,'uint8'),1,bitget(i,samples)); %rotate left    numt = sum(bitget(bitxor(i,j),1:samples));    if numt <= 2      table(i+1) = sum(bitget(i,1:samples));    else      table(i+1) = samples+1;    end  endendmapping.table=table;mapping.samples=samples;mapping.num=newMax;

LBP函数的第五个参数，作者只给出了三种选择，'h'，‘hist’和‘uh’。如果没有设置，程度自动设置为mode='h'。如果想要看一下图像的原始LBP特征，而非直方图，可以随便设一个其它的值。

clear;clc;close allI=imread('test1.bmp');mapping=getmapping(8,'u2');H1=lbp(I,1,8,mapping,'a');%%%如果想以图像形式显示lbp特征，第五个参数随便设一个值，但不能不设。H1=lbp(I,1,8,mapping,1);figure;imshow(H1)H1=lbp(I,1,8,mapping,'h'); %LBP histogram in (8,1) neighborhood                                 %using uniform patternsfigure                                 subplot(2,1,1),stem(H1);H2=lbp(I);subplot(2,1,2),stem(H2);