我有一个图像，我想在上面计算圆形社区内的局部直方图。社区的大小由radius给出。尽管下面的代码可以完成这项工作，但它的计算成本很高。我运行分析器，我访问圆形邻域内像素的方式已经很昂贵了。
是否有任何基于矢量化的改进/优化？或者例如，将邻域存储为列？我在这篇文章中发现了一个类似的问题，建议的解决方案非常符合以下代码的精神，但是该解决方案仍然不适合我的情况。任何想法都非常受欢迎：-）想象一下，目前图像是二进制的，但该方法也应该理想地适用于灰度图像：-）
[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
按照我在回答类似问题时描述的技术，您可以尝试执行以下操作：
- 计算特定体素的索引偏移，从而到达半径内的所有相邻要素
- 确定哪些体素的所有相邻体素距离边的半径至少
- 计算所有这些体素的相邻体素
- 为每个社区生成直方图
对此进行矢量化并不难，但请注意
- 当社区很大时会很慢
- 它涉及生成一个中间矩阵，该矩阵是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
你是
对的，我认为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时，它将捕获更多像素。它会坚持使用你自己的功能，反正速度更快。
您可以尝试使用相反的逻辑进行处理（如注释中简要说明的那样）
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