5

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

让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功能？如何进一步优化我的代码？

提前致谢。

4

2 回答 2

5

正如 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

或者，正如 Milan 所建议的，您可以通过 卷积应用图像膨胀：

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更快，等等）。

于 2013-09-03T11:49:22.003 回答
1

如果只对r==1感兴趣，看makelut和对应的函数bwloolup。

[编辑]

% Let u(i,j) denote the pixels of the matrix u. Calculate the zero crossing
% pixels. These are the pixels in the grid if u(i-1,j) and u(i+1,j) are of
% opposite signs or u(i,j-1) and u(i,j+1) are of opposite signs.

% First, create a function which will us if a pixel is a zero crossing
% pixel, given its 8 neighbors. 

% uSign = u>0; % Note - 0 is treated as negative here. 

% uSign is 3x3, we are evaluating the center pixel uSign(2,2)
zcFun = @(uSign) (uSign(1,2)~=uSign(3,2)) || (uSign(2,1)~=uSign(2,3));

% Make a look up table which tells us what the output should be for the 2^9
% = 512 possible patterns of 3x3 arrays with 1's and 0's.
lut   = makelut(zcFun,3);

% Test image
im = double(imread('pout.tif'));
% Create positve and negative values
im = im -mean(im(:));

% Apply lookup table
imSign = im>0;
imZC   = bwlookup(imSign,lut);

imshowpair(im, imZC);
于 2013-09-03T13:03:57.307 回答