matlab：圆形邻域内局部直方图的有效计算

Question

我有一张图像，我想在其上计算圆形邻域内的局部直方图。邻域的大小由radius 给出。尽管下面的代码可以完成这项工作，但它的计算成本很高。我运行分析器，我访问圆形邻域内的像素的方式已经很昂贵了。

是否有任何基于矢量化的改进/优化？或者例如，将社区存储为列？ 我在post 中发现了一个类似的问题，并且建议的解决方案非常符合下面代码的精神，但是该解决方案仍然不适合我的情况。任何想法都非常受欢迎:-) 现在想象一下，图像是二进制的，但该方法也应该理想地适用于灰度图像:-)

[rows,cols] = size(img);
hist_img      = zeros(rows, cols, 2);
[XX, YY]      = meshgrid(1:cols, 1:rows);
for rr=1:rows
        for cc=1:cols
            distance      = sqrt( (YY-rr).^2 + (XX-cc).^2  );
            mask_radii = (distance <= radius);
            bwresponses   = img(mask_radii);
            [nelems, ~]   = histc(double(bwresponses),0:255);
            % do some processing over the histogram
            ...
        end
end

编辑 1 鉴于收到的反馈，我尝试更新解决方案。但是，它还不正确

radius = sqrt(2.0);
disk   = diskfilter(radius);
fun    = @(x) histc( x(disk>0), min(x(:)):max(x(:)) ); 
output = im2col(im, size(disk), fun);

function disk = diskfilter(radius)
    height  = 2*ceil(radius)+1;
    width   = 2*ceil(radius)+1;
    [XX,YY] = meshgrid(1:width,1:height);
    dist    = sqrt((XX-ceil(width/2)).^2+(YY-ceil(height/2)).^2);
    circfilter = (dist <= radius);
end

Answer 1

按照我在answer to a similar question 中描述的技术，您可以尝试执行以下操作：

1. 计算特定体素的索引偏移量，让您到达半径内的所有邻居
2. 确定哪些体素的所有邻居都至少距离边缘半径
3. 计算所有这些体素的邻居
4. 为每个社区生成直方图

矢量化并不难，但请注意

1. 社区很大时会很慢
2. 它涉及生成一个 NxM 的中间矩阵（N = 图像中的体素，M = 邻域中的体素），它可能会变得非常大

代码如下：

% generate histograms for neighborhood within radius r
A = rand(200,200,200);
radius = 2.5;
tic
sz=size(A);
[xx yy zz] = meshgrid(1:sz(2), 1:sz(1), 1:sz(3));
center = round(sz/2);
centerPoints = find((xx - center(1)).^2 + (yy - center(2)).^2 + (zz - center(3)).^2 < radius.^2);
centerIndex = sub2ind(sz, center(1), center(2), center(3));

% limit to just the points that are "far enough on the inside":
inside = find(xx > radius+1 & xx < sz(2) - radius & ...
    yy > radius + 1 & yy < sz(1) - radius & ...
    zz > radius + 1 & zz < sz(3) - radius);

offsets = centerPoints - centerIndex;
allPoints = 1:prod(sz);
insidePoints = allPoints(inside);
indices = bsxfun(@plus, offsets, insidePoints);

hh = histc(A(indices), 0:0.1:1);  % <<<< modify to give you the histogram you want
toc

相同代码的 2D 版本（这可能就是您所需要的，而且速度要快得多）：

% generate histograms for neighborhood within radius r
A = rand(200,200);
radius = 2.5;
tic
sz=size(A);
[xx yy] = meshgrid(1:sz(2), 1:sz(1));
center = round(sz/2);
centerPoints = find((xx - center(1)).^2 + (yy - center(2)).^2  < radius.^2);
centerIndex = sub2ind(sz, center(1), center(2));

% limit to just the points that are "far enough on the inside":
inside = find(xx > radius+1 & xx < sz(2) - radius & ...
    yy > radius + 1 & yy < sz(1) - radius);

offsets = centerPoints - centerIndex;
allPoints = 1:prod(sz);
insidePoints = allPoints(inside);
indices = bsxfun(@plus, offsets, insidePoints);

hh = histc(A(indices), 0:0.1:1);  % <<<< modify to give you the histogram you want
toc

Answer 2

你是对的，我不认为colfilt 可以使用，因为你没有应用过滤器。您必须检查正确性，但这是我使用im2col 和您的diskfilter 函数的尝试（我确实删除了到double 的转换，所以它现在输出逻辑）：

function circhist

% Example data
im = randi(256,20)-1;

% Ranges - I do this globally for the whole image rather than for each neighborhood
mini = min(im(:));
maxi = max(im(:));
edges = linspace(mini,maxi,20);

% Disk filter
radius = sqrt(2.0);
disk = diskfilter(radius); % Returns logical matrix

% Pad array with -1
im_pad = padarray(im, (size(disk)-1)/2, -1);

% Convert sliding neighborhoods to columns
B = im2col(im_pad, size(disk), 'sliding');

% Get elements from each column that correspond to disk (logical indexing)
C = B(disk(:), :);

% Apply histogram across columns to count number of elements
out = histc(C, edges)

% Display output
figure
imagesc(out)
h = colorbar;
ylabel(h,'Counts');
xlabel('Neighborhood #')
ylabel('Bins')
axis xy

function disk = diskfilter(radius)
height  = 2*ceil(radius)+1;
width   = 2*ceil(radius)+1;
[XX,YY] = meshgrid(1:width,1:height);
dist    = sqrt((XX-ceil(width/2)).^2+(YY-ceil(height/2)).^2);
disk = (dist <= radius);

如果您想根据每个邻域设置范围 (edges)，那么如果您想构建一个大矩阵（然后是该矩阵的行），则需要确保该向量始终具有相同的长度矩阵不会相互对应）。

您应该注意fspecial 返回的磁盘形状不像您使用的那样圆。它旨在用于平滑/平均滤波器，因此边缘是模糊的（抗锯齿）。因此，当您使用~=0 时，它会抓取更多像素。它会坚持使用您自己的功能，反正速度更快。

Answer 3

您可以尝试使用相反的逻辑进行处理（如评论中简要说明的那样）

hist = zeros(W+2*R, H+2*R, Q);
for i = 1:R+1;
  for j = 1:R+1;  
      if ((i-R-1)^2+(j-R-1)^2 < R*R)
         for q = 0:1:Q-1;
             hist(i:i+W-1,j:j+H-1,q+1) += (image == q);
         end
      end
   end
end