图像特征提取——韦伯局部描述符（WLD）

一、原理及概述

　　韦伯局部描述符（WLD）是一种鲁棒性好、简单高效的局部特征描述符。WLD由两个部分组成：差分激励和梯度方向。

　　其具体算法是对于给定的一幅图像，通过对每个像素进行这两个分量的计算来提取其差分激励图像和梯度方向图像，并将这两部分组合成一个二维直方图，进而转化为一个一维直方图，最后利用一维直方图进行分类。^[1]

　　差分激励：差分激励时某像素领域值和减去8倍中心像素值的结果与该像素点的比值。显然，这反映了当前像素点与领域像素点之间的差异。

　　在对差分激励有了一个概念性的认识以后，我们来看数学表述：

　　综合上述公式，我们整合一下：

　　通过计算式子，我们知道如果G_arctan>0，则表明领域像素的亮度大于当前像素；反之，则表明当前像素在该区域内的亮度较低。

　　在计算过程中，为了避免出现输入值过大或过小的情况从而造成不可控的局面，可以通过一个函数将其映射到一个合理范围之内。常用的是sigmoid函数、arctan函数以及对数函数。sigmoid函数如下：

　　由于考虑到sigmoid函数的计算量较大，所以往往可以选用arctan函数。

　　梯度方向：通过计算局部窗内水平方向和垂直方向像素插值之比的反正切变换得到^[1]。梯度方向反映的是局部窗内灰度变化的空间分布信息。

　　以下为梯度方向的数学表述，注意，x_i在下面一张图里。

　　为了更加有效地区分灰度窗内灰度变化的分布情况，，我们需要将其扩展到[0,2pi)。设映射f：，有：

　　最后，需要对其量化，计算过程如下：

　　直方图：在图像中利用直方图表示灰度分布。在此处，直方图是用于局部描述算子的变量之一。

　　我们在得到了差分激励与梯度方向以后，可以得到两幅新的图像：差分激励图像和梯度方向图像，从而可以得到二维WLD直方图。其中，N为图像的维度，T为方向域的分段个数，每一小块都对应一个差分激励和梯度方向。

　　为便于分类，，需要将二维直方图转化为一维直方图。整个过程如下图所示，其实我的理解就是，一行一行的排列，（就像在C里，二维数组其实是个一维数组一样，只是他是一行一行的读取罢了）

　　综合：了解了这几个基本信息以后，我们来看一张图，好详细理解WLD算子：

　　对于差分激励来讲，v₀₀，v₀₁，我们知道他是可以看作上图中的滤波算子，同理，v₁₀，v₁₁也是如此。我们通过分别计算完差分激励与梯度方向后，得到有关二者的二维直方图，最后将二维直方图一维化，即为WLD算子的全部操作。

二、代码及其结果

　　以下为WLD算子的代码：

image=imread('lena.jpg');
image=rgb2gray(image);
d_image=double(image);
[xsize ysize] = size(d_image);
EPSILON=0.0000001;
numNeighbors=8; d_differential_excitation = zeros(xsize-1,ysize-1);
d_gradient_orientation = zeros(xsize-1,ysize-1);%计算WLD
for x = 2:xsize-1for y = 2:ysize-1% 计算差励v00=d_image(x-1,y-1)+d_image(x,y-1)+d_image(x+1,y-1)+d_image(x-1,y)+d_image(x+1,y)+d_image(x-1,y+1)+d_image(x,y+1)+d_image(x+1,y+1)-8*d_image(x,y);v01=d_image(x,y);if ( v01 ~= 0 )d_differential_excitation(x,y)=atan(v00/v01);elsed_differential_excitation(x,y)=0;end% 计算梯度方向N2=d_image(x,y-1);         N6=d_image(x,y+1);N4=d_image(x+1,y);         N0=d_image(x-1,y);if ( abs(N0-N4) < EPSILON)		d_gradient_orientation(x,y) = 0;elsev10=N6-N2;		v11=N0-N4;d_gradient_orientation(x,y) = atan(v11/v10);% 映射至0-2pi之间%角度转化为数值 && 梯度方向的量化d_gradient_orientation(x,y)=d_gradient_orientation(x,y)*180/pi;if     (v11 >  EPSILON && v10 >  EPSILON)		d_gradient_orientation(x,y)= d_gradient_orientation(x,y)+ 180;   else if (v11 < -EPSILON && v10 >  EPSILON)      d_gradient_orientation(x,y)= d_gradient_orientation(x,y)+ 180;  else if (v11 < -EPSILON && v10 < -EPSILON)      d_gradient_orientation(x,y)= d_gradient_orientation(x,y); else if (v11 >  EPSILON && v10 < -EPSILON)		d_gradient_orientation(x,y)= d_gradient_orientation(x,y)+ 360; endendendendendend
end
% 对梯度方向值与差励进行量化
% 将差励6等分,梯度8等分
%reshape函数重新调整矩阵的行数、列数、维数。下为将所有元素展现为一列
d_differential_excitation_one = reshape(d_differential_excitation,[size(d_differential_excitation,1)*size(d_differential_excitation,2),1]);
d_gradient_orientation_one = reshape(d_gradient_orientation,[size(d_gradient_orientation,1)*size(d_gradient_orientation,2),1]);
% figure;
%直方图，指定柱状数量为6
H_dif = histogram(d_differential_excitation_one,6);
dif_edge=H_dif.BinEdges;
% figure;
H_gra = histogram(d_gradient_orientation_one,8);
gra_edge=H_gra.BinEdges;
% count 两个矩阵的像素点的直方图归到6*8矩阵 
%循环实现计数功能，即每一节直方图的出现的像素点的个数
histNumbers = zeros(6,8);
for i=1:length(d_differential_excitation_one)temp1 = d_differential_excitation_one(i);temp2 = d_gradient_orientation_one(i);for k = 2:length(dif_edge)if temp1 <= dif_edge(k)g1 = k-1;break;endendfor k = 2:length(gra_edge)if temp2 <= gra_edge(k)g2 = k-1;break;endendhistNumbers(g1,g2) = histNumbers(g1,g2)+1;
end%figure(111);%bar3(histNumbers');
mid_dif=[];%二维直方图的完善
for i=1:length(dif_edge)-1temp = (dif_edge(i)+dif_edge(i+1))/2;temp = round(temp*100)/100;mid_dif = [mid_dif,temp];
end
set(gca, 'XTickLabel',{num2str(mid_dif(1)),num2str(mid_dif(2)),num2str(mid_dif(3)),num2str(mid_dif(4)),num2str(mid_dif(5)),num2str(mid_dif(6))});
mid_gra=[];
for i=1:length(gra_edge)-1temp = (gra_edge(i)+gra_edge(i+1))/2;temp = round(temp*100)/100;temp = temp*2*pi/360;mid_gra = [mid_gra,temp];
endbarhang=[];
for i=1:size(histNumbers,1)barhang=[barhang,histNumbers(i,:)];
end
% bar(barhang);
% % 绘图
figure(123);
subplot(2,2,1),imshow(d_differential_excitation,[]); title('差励')
image_DE = (255* (d_differential_excitation - (-pi/2))/((pi/2)-(-pi/2)));
subplot(2,2,2),imshow(mat2gray(uint8(image_DE)));title('差励归到[0,255]')
subplot(2,2,3),imshow(d_gradient_orientation,[]); title('梯度方向')
image_DO = (255* (d_gradient_orientation - (-pi/2))/((pi/2)-(-pi/2)));
subplot(2,2,4),imshow(mat2gray(uint8(image_DO)));title('梯度方向归到[0,255]')

　　结果如下：

　　差分激励和梯度方向的两幅图像分别显示如下：