LBP代码详细注释

素颜马尾好姑娘i 2022-06-05 10:20 321阅读 0赞

%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};% % %% 把第一个参数赋值给image
    d_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';
    end
    end
    % % % 如果参数个数大于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
    end
    end
    % % % 如果要参数列表指定了直方图的属性，计算直方图
    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);
    end
    else
    %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);
    end
    end
    
    end
    
     
    
     
    
     
    
    生成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 code
    index = 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
    end
    end
    
    if 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);
    end
    end
    
    if 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
    end
    end
    
    mapping.table=table;
    mapping.samples=samples;
    mapping.num=newMax;
    
     
    
     
    
     
    
    主函数：
    
    clear;clc;close all
    I=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 patterns
    figure 
    subplot(2,1,1),stem(H1);
    H2=lbp(I);
    subplot(2,1,2),stem(H2);