quickest way to calculate neighbourhood of pixels

Question

I have a matrix which consists of positive and negative values. I need to do these things.

Let u(i,j) denote the pixels of the matrix u.

1. 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.
2. Then I need to calculate the narrow band around these zero crossing pixels. The width of the narrow band is (2r+1)X(2r+1) for each pixel. I am taking r=1 so for it I have to actually get the 8 neighborhood pixels of each zero crossing pixels.

I have done this in a program. Please see it below.

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

I am also attaching the isfront and isband functions.

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

The outputs are given below:, as well as the computation time:

%// Computation time

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

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

When I run the code I do get the correct answers but the computation for the tasks is too much to my liking. Is there a better way to do it ? Especially the isband function? How can I optimize my code further ?

Thanks in advance.

Answer 1

If you are only interested in r==1, look at makelut and the corresponding function bwloolup.

[EDIT]

% 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);

Answer 2

As suggested by EitanT, there is at least bwmorph that already does what you want.

If you do no have access to the image processing toolbox, or just insist on doing it yourself:

You can replace the triple-loop in isfront with the vectorized

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;

And you can replace isband by this single loop:

[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

Alternatively, as suggested by Milan, you can apply the image dilation through convolution:

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

which should be even faster.

And you can of course have some other minor optimizations (isband doesn't require phi as an argument, you can pass front as a logical array so that find is faster, etc.).