计算像素邻域的最快方法

Question

我有一个由正值和负值组成的矩阵.我需要做这些事情.

设u(i,j)表示矩阵的像素u.

1. 计算过零像素.这些是在网格如果像素u(i-1,j)和u(i+1,j)具有相反的符号或u(i,j-1)和u(i,j+1)具有相反的符号.
2. 然后我需要计算这些零交叉像素周围的窄带.窄带的宽度(2r+1)X(2r+1)适用于每个像素.我正在考虑r=1它,我必须实际获得每个零交叉像素的8个邻域像素.

我在一个程序中完成了这个.请看下面的内容.

%// calculate the zero crossing pixels  
front = isfront(u);
%// calculate the narrow band of around the zero crossing pixels
band  = isband(u,front,1);

我也附上了isfront和isband功能.

function front = isfront( phi )
%// grab the size of phi
[n, m] = size(phi);

%// create an boolean matrix whose value at each pixel is 0 or 1
%// depending on whether that pixel is a front point or not
front = zeros( size(phi) );

%// A piecewise(Segmentation) linear approximation to the front is contructed by
%// checking each pixels neighbour. Do not check pixels on border.
for i = 2 : n - 1;
  for j = 2 : m - 1;
    if (phi(i-1,j)*phi(i+1,j)<0) || (phi(i,j-1)*phi(i,j+1)<0)
        front(i,j) = 100;
    else
        front(i,j) = 0;
    end
  end
end

function band = isband(phi, front, width)
%// grab size of phi
[m, n] = size(phi);

%// width=r=1;
width = 1;

[x,y] = find(front==100);

%// create an boolean matrix whose value at each pixel is 0 or 1
%// depending on whether that pixel is a band point or not
band = zeros(m, n);

%// for each pixel in phi
for ii = 1:m
  for jj = 1:n
    for k = 1:size(x,1)
        if (ii==x(k)) && (jj==y(k))
            band(ii-1,jj-1) = 100;  band(ii-1,jj) = 100; band(ii-1,jj+1) = 100;
            band(ii  ,jj-1) = 100;  band(ii  ,jj) = 100; band(ii,jj+1) = 100;
            band(ii+1,jj-1) = 100;  band(ii+1,jj) = 100; band(ii+1,jj+1) = 100;
        end
    end
  end
end

输出如下:以及计算时间:

%// Computation time

%// for isfront function
Elapsed time is 0.003413 seconds.

%// for isband function
Elapsed time is 0.026188 seconds.

当我运行代码时,我确实得到了正确的答案,但是任务的计算太多了,不管我喜欢什么.有没有更好的方法呢？特别是isband功能？如何进一步优化我的代码？

提前致谢.

Answer 1

正如EitanT所建议的那样,至少bwmorph已经有了你想做的事情.

如果您无权访问图像处理工具箱,或者只是坚持自己这样做:

您可以isfront使用矢量化替换三重循环

front = zeros(n,m);

zero_crossers = ...
    phi(1:end-2,2:end-1).*phi(3:end,2:end-1) < 0 | ...
    phi(2:end-1,1:end-2).*phi(2:end-1,3:end) < 0;

front([...
                   false(1,m)
    false(n-2,1)  zero_crossers  false(n-2,1)
                   false(1,m)                 ]...
) = 100;

你可以用isband这个循环代替:

[n,m] = size(front);
band = zeros(n,m);
[x,y] = find(front);
for ii = 1:numel(x)
    band(...
        max(x(ii)-width,1) : min(x(ii)+width,n),...
        max(y(ii)-width,1) : min(y(ii)+width,m)) = 1;
end

或者,正如米兰所建议的那样,您可以通过卷积应用图像扩张 :

kernel = ones(2*width+1);    
band = conv2(front, kernel);
band = 100 * (band(width+1:end-width, width+1:end-width) > 0);

哪个应该更快.

当然,您可以进行一些其他的小优化(isband不需要phi作为参数,您可以front作为逻辑数组传递,以便find更快等等).