图像分析：LBP特征解析与代码

矫情吗；* 2023-07-17 06:00 3阅读 0赞

> **LBP**(Local Binary Patterns)，即局部二值模式，是一种描述图像局部空间结构的非参数算子。芬兰Oulu大学的T.Ojala等人于1996年提出这个算子用来分析图像纹理特征，并且描述了它在纹理分类中的强区分能力。LBP算子定义为一种灰度尺度不变的纹理算子，是从局部邻域纹理的普通定义得来的。
> 
> 基本思想是：用中心像素的灰度值作为阈值，与它的邻域相比较得到的二进制码来表述局部纹理特征。
> 
> 在纹理分析方面，LBP算子是最好的纹理描述符之一，它的主要优点有以下几点：
> 
> Ø 通过它的定义可知，LBP算子的灰度尺度不随任何单一变换而变化，因此灰度尺度的鲁棒性好，也就是光照条件下的鲁棒性好；
> 
> Ø 计算速度快。由于它可以通过在小邻域内进行比较操作得到，使得在复杂的实时条件下分析图像成为可能；
> 
> Ø 由于 LBP 算子是一种无参数 (Non-Parametric) 的方法，在应用过程中不需要对它的分布进行预先假设。  
> 
> 
> > *\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_*

*![SouthEast][]*

%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('rice.png');
    %       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.
    error(nargchk(1,5,nargin));
    
    image=varargin{1};
    d_image=double(image);
    
    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
    
    if (nargin == 2) && (length(varargin{2}) == 1)
        error('Input arguments');
    end
    
    if (nargin > 2) && (length(varargin{2}) == 1)
        radius=varargin{2};
        neighbors=varargin{3};
        
        spoints=zeros(neighbors,2);
    
        % Angle step.
        a = 2*pi/neighbors;
        
        for i = 1:neighbors
            spoints(i,1) = -radius*sin((i-1)*a);
            spoints(i,2) = radius*cos((i-1)*a);
        end
        
        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
    
    if (nargin > 1) && (length(varargin{2}) > 1)
        spoints=varargin{2};
        neighbors=size(spoints,1);
        
        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.
    [ysize xsize] = size(image);
    
    
    
    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
    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
    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.
    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
    
    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
    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

%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 1.0 beta
    % Authors: Marko Heikkil?and Timo Ahonen
    % Editor： LinJianmin - Huaqiao University
    % Email:    jianmin1990@outlook.com
    
    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 = bitset(bitshift(i,1,samples),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)) 
        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;
        for j = 1:samples-1
          r = bitset(bitshift(r,1,samples),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,samples),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;

[SouthEast]: /images/20230528/2cb75eafbdbf4522afec6e0503864d1a.png