本地二进制模式的 Matlab 代码

Question

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

[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 代码，您可以 - 理论上，我稍后会讨论 - 沿半径为 R 的圆计算 LBP，采样点为 N，因此 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向左，或等效地乘以 2^(i-1) .然后你用result求和.所以在循环结束时，计算实际上等同于您的 C++ 代码，除了它是在整个图像而不是一个像素上完成的。并且 C++ 代码可以看作是具有半径为 1 和 8 个采样点的相邻圆的 matlab 循环的展开版本。至此，LBP map计算完成，接下来的blocks是对LBP map的额外处理(通过映射表重新映射，可选地计算LBP图像的直方图，而不是LBP图像本身)。

现在，对整个脚本进行一些讨论。这里有一个缺陷隐藏在脚本的末尾。事实上，通过代码，你被限制为 32 个邻居，不能再多了，因为最后 LBP 图像被转换到 int32。 .缺陷是变量 result被分配为双矩阵而不是整数矩阵，所以我真的希望更新时没有近似问题 result稍后在转换为整数时，导致 LBP 中的位发生变化。通常不应该有，因为至少有 52 个精度位(根据 IEEE 754 规范的 wikipedia)。我认为这里有风险……相反，我不知道长固定大小、高效位向量的 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 is better as 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]; .