LBP紋理特征 · 圖像超分辨率技術

局部二進制模式（Local binary patterns，LBP）最早是作為一種有有效的紋理描述算子提出的，由于其對圖像局部紋理特征的卓越描繪能力而獲得了非常廣泛的應用。LBP特征具有很強的分類能力（Highly Discriminative）、較高的計算效率并且對于單調的灰度變化具有不變性。LBP方法在1994年首先由T. Ojala, M.Pietik?inen, 和 D. Harwood 提出，是一個計算機視覺中用于圖像特征分類的重要方法，后來LBP方法與HOG特征分類器聯合使用，改善了一些數據集上的檢測效果。（1）基本LBP 下圖給出了一個基本的LBP算子，應用LBP算子的過程類似于濾波過程中的模板操作。逐行掃描圖像，對于圖像中的每一個像素點，以該點的灰度作為閾值，對周圍3X3的8鄰域進行二值化，按照一定的順序將二值化的結果組成一個8位二進制數，以此二進制數的值（0~255）作為該點的響應。例如，對下圖中的3X3區域的中心點，以其灰度值68作為閾值，對其8鄰域進行二值化，并且從左上點開始按照順時針方向（具體的順序可以任意，只要統一即可）將二值化的結果組成一個二進制數10001011，即十進制的139，作為中心點的響應。在整個逐行掃描過程結束后，會得到一個LBP響應圖像，這個響應圖像的直方圖被稱為LBP統計直方圖，或LBP直方圖，它常常被作為后續識別工作的特征，因此也被稱為LBP特征。 LBP的主要思想是以某一點與其鄰域像素的相對灰度作為響應，正是這種相對機制使得LBP算子對于單調的灰度變化具有不變性。人臉圖像常常會受到光照因素的影響而產生灰度變化，但在一個局部區域內，這種變化常常可以被視為是單調的，因此LBP在光照不均的人臉識別應用中也取得了很好的效果。 ![](https://box.kancloud.cn/2016-04-07_5706179104dd4.jpg) （2）圓形鄰域LBP算子基本LBP算子可以被進一步推廣為使用不同大小和形狀的鄰域。采用圓形的鄰域并結合雙線性插值運算可以獲得任意半徑和任意數目的鄰域像素點。如下圖，應用LBP算法的三個鄰域示例所示）進行順時針或逆時針的比較，如果中心像素值比該鄰點大，則將鄰點賦值為1，否則賦值為0，這樣每個點都會獲得一個8位二進制數（通常轉換為十進制數）。然后計算每個cell的直方圖，即每個數字（假定是十進制數）出現的頻率（也就是一個關于每一個像素點是否比鄰域內點大的一個二進制序列進行統計），然后對該直方圖進行歸一化處理。最后將得到的每個cell的統計直方圖進行連接，就得到了整幅圖的LBP紋理特征，然后便可利用SVM或者其他機器學習算法進行分類了。 ![](https://box.kancloud.cn/2016-04-07_5706179124f1c.jpg) （３）ＭＡＴＬＡＢ實現一共有三個m文件，一個是lbp.m, 存放主要的lbp算法，一個是getmapping,用以做算法的輔助函數，一個是cont.m。１.lbp.m ~~~ %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.3 % 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 = roundn((1 - tx) * (1 - ty),-6); w2 = roundn(tx * (1 - ty),-6); w3 = roundn((1 - tx) * ty,-6) ; % w4 = roundn(tx * ty,-6) ; w4 = roundn(1 - w1 - w2 - w3, -6); % 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); N = roundn(N,-4); 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 function x = roundn(x, n) error(nargchk(2, 2, nargin, 'struct')) validateattributes(x, {'single', 'double'}, {}, 'ROUNDN', 'X') validateattributes(n, ... {'numeric'}, {'scalar', 'real', 'integer'}, 'ROUNDN', 'N') if n < 0 p = 10 ^ -n; x = round(p * x) / p; elseif n > 0 p = 10 ^ n; x = p * round(x / p); else x = round(x); end end ~~~ ２.getmapping．ｍ ~~~ %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.2 % Authors: Marko Heikkil?, Timo Ahonen and Xiaopeng Hong % 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. % Modified by Xiaopeng HONG and Guoying ZHAO % Changelog % 0.2 % Solved the compatible issue for the bitshift function in Matlab % 2012 & higher matlab_ver = ver('MATLAB'); matlab_ver = str2double(matlab_ver.Version); if matlab_ver < 8 mapping = getmapping_ver7(samples,mappingtype); else mapping = getmapping_ver8(samples,mappingtype); end end function mapping = getmapping_ver7(samples,mappingtype) disp('For Matlab version 7.x and lower'); 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; end function mapping = getmapping_ver8(samples,mappingtype) disp('For Matlab version 8.0 and higher'); 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 i_bin = dec2bin(i,samples); j_bin = circshift(i_bin',-1)'; %circularly rotate left numt = sum(i_bin~=j_bin); %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_bin = dec2bin(i,samples); for j = 1:samples-1 r = bin2dec(circshift(r_bin',-1*j)'); %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 i_bin = dec2bin(i,samples); j_bin = circshift(i_bin',-1)'; numt = sum(i_bin~=j_bin); 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; end ~~~ ３．ｃｏｎｔ．ｍ ~~~ %C computes the VAR descriptor. % J = CONT(I,R,N,LIMS,MODE) returns either a rotation invariant local % variance (VAR) image or a VAR histogram of the image I. The VAR values % are determined for all pixels having neighborhood defined by the input % arguments. The VAR operator calculates variance on a circumference of % R radius circle. The circumference is discretized into N equally spaced % sample points. Function returns descriptor values in a continuous form or % in a discrete from if the quantization limits are defined in the argument % LIMS. % % Examples % -------- % % im = imread('rice.png'); % c = cont(im,4,16); % d = cont(im,4,16,1:500:2000); % % figure % subplot(121),imshow(c,[]), title('VAR image') % subplot(122),imshow(d,[]), title('Quantized VAR image') function result = cont(varargin) % Version: 0.1.0 % 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; lims=0; mode='i'; end if (nargin > 2) && (length(varargin{2}) == 1) radius=varargin{2}; neighbors=varargin{3}; spoints=zeros(neighbors,2); lims=0; mode='i'; % 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 && ~ischar(varargin{4})) lims=varargin{4}; end if(nargin >= 4 && ischar(varargin{4})) mode=varargin{4}; end if(nargin == 5) mode=varargin{5}; end end if (nargin == 2) && ischar(varargin{2}) mode=varargin{2}; spoints=[-1 -1; -1 0; -1 1; 0 -1; -0 1; 1 -1; 1 0; 1 1]; neighbors=8; lims=0; 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; %Compute the local contrast for i = 1:neighbors y = spoints(i,1)+origy; x = spoints(i,2)+origx; % Calculate floors and ceils for the x and y. fy = floor(y); cy = ceil(y); fx = floor(x); cx = ceil(x); % 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); % Compute the variance using on-line algorithm % ( http://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#On-line_algorithm ). if i == 1 MEAN=zeros(size(N)); DELTA=zeros(size(N)); M2=zeros(size(N)); end DELTA=N-MEAN; MEAN=MEAN+DELTA/i; M2=M2+DELTA.*(N-MEAN); end % Compute the variance matrix. % Optional estimate for variance: % VARIANCE_n=M2/neighbors; result=M2/(neighbors-1); % Quantize if LIMS is given if lims [q r s]=size(result); quant_vector=q_(result(:),lims); result=reshape(quant_vector,q,r,s); if strcmp(mode,'h') % Return histogram result=hist(result, length(lims)-1); end end if strcmp(mode,'h') && ~lims % Return histogram %epoint = round(max(result(:))); result=hist(result(:),0:1:1e4); end end function indx = q_(sig,partition) [nRows, nCols] = size(sig); indx = zeros(nRows, nCols); for i = 1 : length(partition) indx = indx + (sig > partition(i)); end end ~~~ ~~~ %%LBP_test.m Img=imread('lena.jpg'); I=rgb2gray(Img); 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. % show the images figure, imshow(I); title('Input Image'); figure, imshow(I2); title('Result of LBP'); ~~~ ![](https://box.kancloud.cn/2016-04-07_570617913812b.jpg) ![](https://box.kancloud.cn/2016-04-07_5706179154423.jpg) ![](https://box.kancloud.cn/2016-04-07_570617917ce00.jpg)