局部二进制模式的Matlab代码

Question

我正在努力理解这里发现的LBP算法的Matlab实现.我试图找到它如何计算每个像素的二进制文件？它只计算相邻像素大于实际中心像素大小的位置.我想计算每个像素的二进制数,以便使用局部直方图来计算图像的特征.

[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
% the whole process here
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
 size(result)
end

迭代地,它为每个像素的所有8个邻居在结果中添加了一些值.但它是如何与LBP二进制相关联的？它与以下c ++ LBP方法的代码如何相关:

 uchar lbp(const Mat_<uchar> & img, int x, int y)
 {
  // this is pretty much the same what you already got..
  uchar v = 0;
  uchar c = img(y,x);
  v += (img(y-1,x  ) > c) << 0;
  v += (img(y-1,x+1) > c) << 1;
  v += (img(y  ,x+1) > c) << 2;
  v += (img(y+1,x+1) > c) << 3;
  v += (img(y+1,x  ) > c) << 4;
  v += (img(y+1,x-1) > c) << 5;
  v += (img(y  ,x-1) > c) << 6;
  v += (img(y-1,x-1) > c) << 7;
  return v;

}

Answer 1

它是LBP的矢量化实现,非常适合Matlab.

在初始化指令之后,让我们看看主循环,从" for i = 1:neighbors" 行开始.循环非常清楚:它计算一个邻居与中心像素的比较,循环迭代所有邻居.你有这一点,所以现在进入循环深入了解它如何积累所有结果.

事实上,循环的核心过于复杂,因为它考虑了实际圆而不是近似整数圆.因此,指令的主要部分的目的是计算相邻像素的内插强度.这里它与您作为参考的C++代码不同,它只采用整数,1像素宽的半径圆.请记住,使用lbp.m代码,您可以 - 理论上,我将在稍后讨论 - 沿着具有N个采样点的半径R的圆计算LBP,因此C++将对应于半径为1的圆并且具有8个采样点,如果只有没有插值.但是当邻居不适合图像的像素网格时,存在插值,当(abs(x - rx) < 1e-6) && (abs(y - ry) < 1e-6)为假时).

如果(abs(x - rx) < 1e-6) && (abs(y - ry) < 1e-6)为真,则没有插值,因此中心像素和当前邻居之间的所有比较的计算直接存储到D.否则,它计算整个图像上采样邻点处的强度的双线性插值: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);.

最后,转到更新部分:v = 2^(i-1); result = result + v*D;.v相当于移位:对于第i个邻居,你将比较的值i-1移到左边,或者相等地乘以be 2^(i-1).然后你总结result.所以在循环结束时,计算实际上等同于你的C++代码,除了它是在整个图像而不是单个像素上完成的.并且C++代码可以看作是带有半径为1和8个采样点的相邻圆的matlab循环的展开版本.此时,计算LBP图,以下块是LBP图的附加处理(通过映射表重新映射,并且可选地计算LBP图像的直方图而不是LBP图像本身).

现在,关于整个脚本的一点讨论.这里有一个隐藏在脚本末尾的缺陷.事实上,通过代码,你只能被限制为32个邻居,因为最终LBP图像被转换为int32.缺点是变量result被分配为双矩阵而不是整数矩阵,所以我真的希望在更新时result以及稍后在转换为整数时没有近似问题,导致LBP中的位改变.通常不应该有至少52个精度位(根据维基百科的IEEE 754规范).我认为这有风险......相反,我不知道长固定大小的有效位向量的matlab类型.我会使用int64而不是int32,但限制将在64个采样邻居.

编辑

现在,如果你想要在3*3邻域上限制一些局部二进制模式,那么这个Matlab函数对你来说太通用了,最好的办法是展开这个邻域的循环,因此非常接近于C++代码.这是一段代码(我使用按位或代替添加,但它是等效的):

result = uint8(ysize, xsize);
result = (image(1:end-2,2:end-1) > image(2:end-1,2:end-1));                                 % <=> v += (img(y-1,x  ) > c) << 0;
result = result|bitshift((image(1:end-2,3:end) > image(2:end-1,2:end-1)), 1, 'uint8');      % <=> v += (img(y-1,x+1) > c) << 1;
result = result|bitshift((image(2:end-1,3:end) > image(2:end-1,2:end-1)), 2, 'uint8');      % <=> v += (img(y  ,x+1) > c) << 2;
result = result|bitshift((image(3:end,3:end) > image(2:end-1,2:end-1)), 3, 'uint8');        % <=> v += (img(y+1,x+1) > c) << 3;
result = result|bitshift((image(3:end,2:end-1) > image(2:end-1,2:end-1)), 4, 'uint8');      % <=> v += (img(y+1,x  ) > c) << 4;
result = result|bitshift((image(3:end,1:end-2) > image(2:end-1,2:end-1)), 5, 'uint8');      % <=> v += (img(y+1,x-1) > c) << 5;
result = result|bitshift((image(2:end-1,3:end) > image(2:end-1,2:end-1)), 6, 'uint8');      % <=> v += (img(y  ,x-1) > c) << 6;
result = result|bitshift((image(1:end-2,1:end-2) > image(2:end-1,2:end-1)), 7, 'uint8');    % <=> v += (img(y-1,x-1) > c) << 7;

这是使用强大的矢量化将C代码精确转换为Matlab脚本.有了这个,在这个社区改变另一个订单或不同的测试非常简单.我也提到这一点,因为在这个案例的Matlab脚本中有一个错误,第53行有一个错误的标志:neighobrhood更好spoints=[-1 -1; -1 0; -1 1; 0 -1; 0 -1; 1 -1; 1 0; 1 1];而不是spoints=[-1 -1; -1 0; -1 1; 0 -1; -0 1; 1 -1; 1 0; 1 1];.

Answer 2

我做了本地二进制模式的最后一年项目.我看到你指向的代码,但我决定编写自己的代码.

这是我的代码:

function [ LBP ] = LBP( I2)
m=size(I2,1);
n=size(I2,2);
for i=2:m-1
    for j=2:n-1
        J0=I2(i,j);
        I3(i-1,j-1)=I2(i-1,j-1)>J0;
        I3(i-1,j)=I2(i-1,j)>J0;
        I3(i-1,j+1)=I2(i-1,j+1)>J0; 
        I3(i,j+1)=I2(i,j+1)>J0;
        I3(i+1,j+1)=I2(i+1,j+1)>J0; 
        I3(i+1,j)=I2(i+1,j)>J0; 
        I3(i+1,j-1)=I2(i+1,j-1)>J0; 
        I3(i,j-1)=I2(i,j-1)>J0;
        LBP(i,j)=I3(i-1,j-1)*2^7+I3(i-1,j)*2^6+I3(i-1,j+1)*2^5+I3(i,j+1)*2^4+I3(i+1,j+1)*2^3+I3(i+1,j)*2^2+I3(i+1,j-1)*2^1+I3(i,j-1)*2^0;
    end
end
end

I2是您传递的图像,LBP是输出.为Local Binary Pattern写了这个:http://quantgreeks.com/local-binary-pattern-in-matlab/.我知道我可以用更有效的形式编写代码.但是以这种方式编写它会清楚地表明本地二进制模式是如何工作的.