Steve on Image Processing

Blog reader DKS asked recently why values outside the image are assumed to be -Inf when computing dilation. I thought this issue was worth exploring further because it has practical implications for certain computations.

Suppose we start with a simple 4-by-4 matrix:

f = [22 23 15 16; 24 25 14 15; 20 18 17 23; 19 16 15 20]

f =

    22    23    15    16
    24    25    14    15
    20    18    17    23
    19    16    15    20

Now think about computing the erosion with a 3-by-3 flat structuring element. What's the output at the upper left corner? It's the minimum of these 9 values:

?? ?? ??
?? 22 23
?? 24 25

So what do we use for the unknown values outside the image boundary? Suppose we zero pad:

 0  0  0
 0 22 23
 0 24 25

Then the (1,1) output value is 0. In fact, if the image pixels are nonnegative, which is common, there will be a zero-valued one-pixel-wide border all the way around the edge of the output image. This is called a boundary artifact and is undesirable. We could avoid this problem by simply excluding external values in our computation of the minimum. A mathematical equivalent for erosion is to assume external values all equal some constant that is guaranteed to be greater than or equal to all image pixels - like Inf:

Inf Inf Inf
Inf  22  23
Inf  24  25

To avoid boundary artifacts when performing dilation, pad with -Inf:

-Inf -Inf -Inf
-Inf   22   23
-Inf   24   25

Most of the time you don't need to explicitly use the Inf values. One type of exception occurs when the structuring element doesn't include the center element. Then sometimes the Inf values can appear in the output.

imdilate(f, [1 0 0])

ans =

    23    15    16  -Inf
    25    14    15  -Inf
    18    17    23  -Inf
    16    15    20  -Inf

Similarly, most of the time you don't actually need to explicitly pad the image. The exception is when you are computing a sequence of dilations (or erosions). In the Image Processing Toolbox this frequently occurs because imdilate and imerode exploit ''structuring element decomposition.'' That is, a large structuring element is decomposed into two or more smaller structuring elements that are mathematically equivalent to the original.

Here's a contrived example to demonstrate why you need to pad explicitly to make the mathematical equivalence work out.

Structuring element 1 translates an image to the left by 2 pixels.

se1 = [1 0 0 0 0];

Structuring element 2 translates an image to the right by 1 pixel.

se2 = [0 0 1];

You'd expect the composition of dilation with these two structuring elements to be equivalent to a single dilation with a structuring element that translates an image left by 1 pixel:

se3 = [1 0 0];

Let's try that sequence with no padding:

f1 = imdilate(f, se1)

f1 =

    15    16  -Inf  -Inf
    14    15  -Inf  -Inf
    17    23  -Inf  -Inf
    15    20  -Inf  -Inf

f2 = imdilate(f1, se2)

f2 =

  -Inf    15    16  -Inf
  -Inf    14    15  -Inf
  -Inf    17    23  -Inf
  -Inf    15    20  -Inf

Now dilate the original image with se3 and compare:

f3 = imdilate(f, se3)

f3 =

    23    15    16  -Inf
    25    14    15  -Inf
    18    17    23  -Inf
    16    15    20  -Inf

They aren't the same. To make the results equivalent, you have to pad the original image, perform the dilation sequence, and the crop the result.

fp = padarray(f, [2 2], -Inf)

fp =

  -Inf  -Inf  -Inf  -Inf  -Inf  -Inf  -Inf  -Inf
  -Inf  -Inf  -Inf  -Inf  -Inf  -Inf  -Inf  -Inf
  -Inf  -Inf    22    23    15    16  -Inf  -Inf
  -Inf  -Inf    24    25    14    15  -Inf  -Inf
  -Inf  -Inf    20    18    17    23  -Inf  -Inf
  -Inf  -Inf    19    16    15    20  -Inf  -Inf
  -Inf  -Inf  -Inf  -Inf  -Inf  -Inf  -Inf  -Inf
  -Inf  -Inf  -Inf  -Inf  -Inf  -Inf  -Inf  -Inf

fp1 = imdilate(fp, se1)

fp1 =

  -Inf  -Inf  -Inf  -Inf  -Inf  -Inf  -Inf  -Inf
  -Inf  -Inf  -Inf  -Inf  -Inf  -Inf  -Inf  -Inf
    22    23    15    16  -Inf  -Inf  -Inf  -Inf
    24    25    14    15  -Inf  -Inf  -Inf  -Inf
    20    18    17    23  -Inf  -Inf  -Inf  -Inf
    19    16    15    20  -Inf  -Inf  -Inf  -Inf
  -Inf  -Inf  -Inf  -Inf  -Inf  -Inf  -Inf  -Inf
  -Inf  -Inf  -Inf  -Inf  -Inf  -Inf  -Inf  -Inf

fp2 = imdilate(fp1, se2)

fp2 =

  -Inf  -Inf  -Inf  -Inf  -Inf  -Inf  -Inf  -Inf
  -Inf  -Inf  -Inf  -Inf  -Inf  -Inf  -Inf  -Inf
  -Inf    22    23    15    16  -Inf  -Inf  -Inf
  -Inf    24    25    14    15  -Inf  -Inf  -Inf
  -Inf    20    18    17    23  -Inf  -Inf  -Inf
  -Inf    19    16    15    20  -Inf  -Inf  -Inf
  -Inf  -Inf  -Inf  -Inf  -Inf  -Inf  -Inf  -Inf
  -Inf  -Inf  -Inf  -Inf  -Inf  -Inf  -Inf  -Inf

fp2_cropped = fp2(3:end-2, 3:end-2)

fp2_cropped =

    23    15    16  -Inf
    25    14    15  -Inf
    18    17    23  -Inf
    16    15    20  -Inf

isequal(f3, fp2_cropped)

ans =

     1

Whenever imdilate or imerode is called with a decomposed structuring element, the functions compute the minimum padding necessary to avoid boundary artifacts and pad with either -Inf or Inf, respectively.

It's a classic speed vs. memory tradeoff. Exploiting structuring element decomposition is faster, but storing the intermediate padded arrays takes more memory.

Get the MATLAB code

Published with MATLAB® 7.4