A fascinating question came to the Image Processing Toolbox development team recently from Brett Shoelson, a MathWorks application engineer, prolific File Exchange contributor, and former MATLAB Central blogger. Brett was trying to solve a problem for a customer. Part of the problem came down to this question:

Suppose you have a label matrix that represents a number of different regions. Are any of the regions disconnected?

Let me elaborate. A label matrix looks something like this:

L = [ ...
    1 1 1 2 5 6
    1 1 2 2 5 5
    3 4 4 4 4 4
    3 3 7 7 7 8
    3 3 7 7 8 8
    9 9 9 9 9 4]

L =

     1     1     1     2     5     6
     1     1     2     2     5     5
     3     4     4     4     4     4
     3     3     7     7     7     8
     3     3     7     7     8     8
     9     9     9     9     9     4

For a positive integer k, the set of elements of L equal to k is the kth region. So, for example, the 5th region is:

L == 5

ans =

     0     0     0     0     1     0
     0     0     0     0     1     1
     0     0     0     0     0     0
     0     0     0     0     0     0
     0     0     0     0     0     0
     0     0     0     0     0     0

Or, displayed as an image:

imshow(L == 5,'InitialMagnification','fit')

So my sample label matrix here represents 9 different regions. One of them, region 4 is disconnected, meaning it consists of more than one connected set (or connected component) or pixels.

L == 4

ans =

     0     0     0     0     0     0
     0     0     0     0     0     0
     0     1     1     1     1     1
     0     0     0     0     0     0
     0     0     0     0     0     0
     0     0     0     0     0     1

imshow(L == 4,'InitialMagnification','fit')
title('Region 4')

How can we programmatically determine whether any of the regions is disconnected? A brute force method is to compute, one at a time, the binary image associated with each region and then call bwconncomp to determine the number of connected components associated with it.

for k = 1:9
    mask_k = (L == k);
    cc = bwconncomp(mask_k);
    if cc.NumObjects > 1
        fprintf('Region %d is disconnected.\n',k)
    end
end

Region 4 is disconnected.

If a brute force method gives you the desired answer and does it fast enough, then sometimes you just declare success and move on to solve the next problem of the day. Was this method fast enough? "No," said Brett. "My customer's images are about 25,000-by-25,000, and this method is way too slow."

"Hmm," I replied, trying to sound thoughtful while really just playing for time.

"OK, ask your customer this question: Is it known whether any of the regions in question contain any holes?"

Soon, the answer came back: "No, the labeled regions do not contain holes."

That was the answer I had hoped for. I started thinking about how to compute something called the Euler number for every labeled region, more or less simultaneously.

No, not that Euler number. I was thinking of the Euler number for a binary image, which is the number of objects in the image (or the number of connected components) minus the number of holes. (I don't know why this property of a binary image is called Euler number. If you know, please let us know in the comments.) The Image Processing Toolbox function bweuler computes this number for a binary image.

For a topological property of an image, the Euler number is unusual because it can be computed by one sweep of a purely local, 2-by-2 neighborhood operation. I first saw the method in the book Robot Vision, Berthold K. P. Horn, MIT Press, 1989. In this method, you count the number of occurrences of this 2-by-2 neighborhood pattern (called an upstream convexity):

0 0
0 1

and then subtract the number of occurrences of this 2-by-2 neighborhood pattern (called an upstream concavity):

0 1
1 1

(Note: this variation of the method works for 4-connected objects, not for 8-connected objects.)

But Brett has a label matrix representing many different regions, not a single binary image mask representing one region. How can we apply this technique?

To answer the question, let's look at the 2-by-2 neighborhood to the northwest of element (3,2):

L32 = L(2:3,1:2)

L32 =

     1     1
     3     4

I'll take the lower-right pixel as my reference point for the neighborhood. I can look at just the binary mask formed by comparing L(3,2) to its three neighbors to the north and west.

L(3,2) == L32

ans =

     0     0
     0     1

So, for region 4, we have an upstream convexity at this location. We could repeat this calculation in one loop over all the pixel locations, looking for (and counting) upstream convexities and upstream concavities associated with the label at the reference point for each 2-by-2 neighborhood.

I want to show you how to do this just with matrix operations and indexing. No loop is required. The technique is fairly general and can be applied to other local neighborhood operations.

Consider these two submatrices:

L1 = L(1:end-1,1:end-1)

L1 =

     1     1     1     2     5
     1     1     2     2     5
     3     4     4     4     4
     3     3     7     7     7
     3     3     7     7     8

and

L2 = L(2:end,2:end)

L2 =

     1     2     2     5     5
     4     4     4     4     4
     3     7     7     7     8
     3     7     7     8     8
     9     9     9     9     4

Note that the elements in L1 are the northwest neighbors of the corresponding elements in L2. This basic idea gives us a general approach for doing certain kinds of neigbhorhood operations without explicit per-pixel loops.

Here's a short function I wrote that uses this technique to compute 4-connected Euler numbers for a label matrix:

type eulerNumbers4

function e = eulerNumbers4(L)
% e = eulerNumbers4(L)
% 
% For a label matrix L containing nonnegative integers, returns a vector e
% such that e(k) is the 4-connected Euler number of the binary image (L ==
% k), for k = 1:max(L(:)).

Lp = padarray(L,[1 1],0,'pre');
num_regions = max(Lp(:));

I_NW = Lp(1:end-1,1:end-1);
I_N = Lp(1:end-1,2:end);
I_W = Lp(2:end,1:end-1);

is_upstream_convexity = (L > 0) & (L ~= I_N);
is_upstream_convexity = is_upstream_convexity & (L ~= I_NW);
is_upstream_convexity = is_upstream_convexity & (L ~= I_W);

is_upstream_concavity = (L > 0) & (L ~= I_NW);
is_upstream_concavity = is_upstream_concavity & (L == I_N);
is_upstream_concavity = is_upstream_concavity & (L == I_W);

upstream_convexity_labels = L(is_upstream_convexity);
upstream_concavity_labels = L(is_upstream_concavity);

total_upstream_convexities = accumarray(upstream_convexity_labels, ...
    1,[num_regions 1]);
total_upstream_concavities = accumarray(upstream_concavity_labels, ...
    1,[num_regions 1]);

e = total_upstream_convexities - total_upstream_concavities;

Under the assumption that the regions do not contain holes, the Euler number for each region is exactly the number of 4-connected components belong to that region.

eulerNumbers4(L)

ans =

     1
     1
     1
     2
     1
     1
     1
     1
     1

So, we can finally answer the question: Are any of the labeled regions disconnected?

any(eulerNumbers4(L) > 1)

ans =

     1

Yes.

True confessions, though. After I sent this version to Brett, and it didn't work for his problem, we realized that he needed the answer for 8-connected regions, not 4-connected.

It turns out that there is a local neighborhood-based Euler number algorithm for 8-connected regions, but it is significantly more complex to implement in this form. I did write it up, but I thought it would be too much uninteresting code to show here. If you are really interested, the method is described in Pratt, Digital Image Processing, 2nd edition, John Wiley & Sons, 1991, pp. 632-634. Pratt, in turn, credits the paper S. B. Gray, "Local Properties of Binary Images in Two Dimensions," IEEE Transactions on Computers, C-20, 5, May 1971, pp. 551-561.

The function bweuler is based on this algorithm, plus the binary image neighborhood lookup table approach that I have described here previously.

Have you used the Euler number measurement in your own work? Let us know in the comments.

Get the MATLAB code

Published with MATLAB® R2014a