Steve on Image Processing

A friend from my grad school days (back in the previous millenium) is an electrical engineering professor. Students in his class recently asked him why the conv function is not implemented using FFTs.

I'm not on the team responsible for conv, but I wrote back with my thoughts, and I thought I would share them here as well.

Let's review the basics. Using the typical convolution formula to compute the one-dimensional convolution of a P-element sequence A with Q-element sequence B has a computational complexity of . However, the discrete Fourier transform (DFT) can be used to implement convolution as follows:

1. Compute the L-point DFT of A, where .

2. Compute the L-point DFT of B.

3. Multiply the two DFTs.

4. Compute the inverse DFT to get the convolution.

Here's a simple MATLAB function for computing convolution using the Fast Fourier Transform (FFT), which is simply a fast algorithm for computing the DFT.

type conv_fft

function c = conv_fft(a, b)

P = numel(a);
Q = numel(b);
L = P + Q - 1;
K = 2^nextpow2(L);

c = ifft(fft(a, K) .* fft(b, K));
c = c(1:L);

Note that the code uses the next power-of-two greater than or equal to L, although this is not strictly necessary. The fft function operates in time whether or not L is a power of two.

The overall computational complexity of these steps is . For P and Q sufficiently large, then, using the DFT to implement convolution is a computational win.

So why don't we do it?

There are several technical factors. Let's look at speed, exact computation, and memory.

Speed

One factor is that DFT-based computation is not always faster. Let's do an experiment where we compute the convolution of a 1000-element sequence with another sequence of varying length. (Get the timeit function from the MATLAB Central File Exchange.)

x = rand(1, 1000);
nn = 25:25:1000;

t_normal = zeros(size(nn));
t_fft = zeros(size(nn));

for k = 1:numel(nn)
    n = nn(k);
    y = rand(1, n);
    t_normal(k) = timeit(@() conv(x, y));
    t_fft(k) = timeit(@() conv_fft(x, y));
end

plot(nn, t_normal, nn, t_fft)
legend({'Normal computation', 'FFT-based computation'})

For sequences y shorter than a certain length, called the cross-over point, it's quicker to use the normal computation.

Exact computation

A second consideration is whether the computation is subject to floating-point round-off errors and to what degree. There are applications, for example, where the convolution of integer-valued sequences is computed. For such applications, a user would reasonably expect the output sequence to be integer-valued as well.

To illustrate, here's a simple function that computes n-th order binomial coefficients using convolution:

type binom

function c = binom(n)

c = 1;
for k = 1:n
    c = conv(c, [1 1]);
end

binom(7)

ans =

     1     7    21    35    35    21     7     1

I wrote a variation called binom_fft that is the same as binom except that it calls conv_fft.

format long
binom_fft(7)

ans =

  Columns 1 through 4

   1.000000000000000   6.999999999999998  21.000000000000000  35.000000000000000

  Columns 5 through 8

  35.000000000000000  21.000000000000000   7.000000000000000   0.999999999999996

Whoops! What's going on? The answer is that the FFT-based implementation of convolution is subject to floating-point round-off error.

I imagine that most MATLAB users would consider the output of binom_fft to be wrong.

Memory

The last technical consideration I want to mention is memory. Because of the padding and complex arithmetic involved in the FFT computations, the FFT-based convolution implementation requires a lot more memory than the normal computation. This may not often be a problem for one-dimensional computations, but it can be a big deal for multidimensional computations.

Final thoughts

The technical considerations listed above can all be solved in principle. The implementation could switch to using the normal method for short or integer-valued sequences, for example. And there are FFT-based techniques such as overlap-and-add to reduce the memory load.

But the problem can get quite complicated. Testing floating-point values to see if they are integers, for example, can be slow. Also, I suspect (but have not checked) that a multithreaded implemention of the normal computation will take advantage of multiple cores better than a multithreaded FFT-based method. That would change the cross-over point between the two methods. And the exact cross-over point will vary from computer to computer, making it likely that our implementation would be somewhat slower for some people for some problems.

Overall, my inclination would be to provide FFT-based convolution as a separate function rather than reimplementing conv. And that's what we've done. See the function fftfilt in the Signal Processing Toolbox.

Do you disagree with this approach? Post your comment.

Get the MATLAB code

Published with MATLAB® 7.9

The R2009b release of MATLAB contains a new Tiff class. The primary purpose of this class is to provide lower-level access to creating and modifying image data and metadata in an existing TIFF file.

Several customers have asked for the ability to embed an ICC profile into a TIFF file. The Image Processing Toolbox function iccread can read a profile embedded in a TIFF file, but iccwrite can only write a stand-alone profile file.

Here is code using the new Tiff class to embed a profile in an existing TIFF file. Let's start by rewriting one of the Image Processing Toolbox sample PNG image files as a TIFF file.

rgb = imread('peppers.png');
imwrite(rgb, 'peppers.tif');
s = dir('peppers.tif')

s = 

       name: 'peppers.tif'
       date: '20-Oct-2009 09:14:29'
      bytes: 593880
      isdir: 0
    datenum: 7.3407e+005

Now let's embed the sample profile sRGB.icm into the TIFF file we just made.

Step 1. Read in the raw bytes of the profile file.

fid = fopen('sRGB.icm');
raw_profile_bytes = fread(fid, Inf, 'uint8=>uint8');
fclose(fid);

Step 2. Initialize a Tiff object using 'r+' mode (read and modify).

tif = Tiff('peppers.tif', 'r+')

tif = 

                  TIFF File: 'peppers.tif'
                       Mode: 'r+'
    Current Image Directory: 1
           Number Of Strips: 77
                SubFileType: Tiff.SubFileType.Default
                Photometric: Tiff.Photometric.RGB
                ImageLength: 384
                 ImageWidth: 512
               RowsPerStrip: 5
              BitsPerSample: 8
                Compression: Tiff.Compression.PackBits
               SampleFormat: Tiff.SampleFormat.UInt
            SamplesPerPixel: 3
        PlanarConfiguration: Tiff.PlanarConfiguration.Chunky
                Orientation: Tiff.Orientation.TopLeft

Step 3. Embed the profile bytes as a TIFF tag.

tif.setTag('ICCProfile', raw_profile_bytes);

Step 4. Tell the Tiff object to update the image metadata in the file.

tif.rewriteDirectory();

Step 5. Close the Tiff object.

tif.close();

Now the TIFF file contains the profile. Notice the file size has changed.

s = dir('peppers.tif')

s = 

       name: 'peppers.tif'
       date: '20-Oct-2009 09:14:30'
      bytes: 597828
      isdir: 0
    datenum: 7.3407e+005

And we can read in the profile using iccread.

p = iccread('peppers.tif')

p = 

               Header: [1x1 struct]
             TagTable: {17x3 cell}
            Copyright: 'Copyright (c) 1999 Hewlett-Packard Company'
          Description: [1x1 struct]
      MediaWhitePoint: [0.9505 1 1.0891]
      MediaBlackPoint: [0 0 0]
        DeviceMfgDesc: [1x1 struct]
      DeviceModelDesc: [1x1 struct]
      ViewingCondDesc: [1x1 struct]
    ViewingConditions: [1x1 struct]
            Luminance: [76.0365 80 87.1246]
          Measurement: [1x1 struct]
           Technology: 'Cathode Ray Tube Display'
               MatTRC: [1x1 struct]
          PrivateTags: {}
             Filename: 'peppers.tif'

It would be logical to enhance iccwrite to make this easier for you. We'll look into that, although I'm not sure when that might happen.

Get the MATLAB code

Published with MATLAB® 7.9

We've been working for a while now to make Image Processing Toolbox functions run faster. The R2009b release notes mention several performance improvements. We've gotten some feedback, though, that our release notes are pretty vague about the improvements. I can't argue with that impression. We tend to be vague because performance optimization is a very complex topic, and it can be quite difficult to characterize performance changes in a way that is brief, understandable, and accurate for every user's own hardware and data.

But I've decided to start posting more detailed information here about the performance improvements. I have more flexibility (and room!) here than we have with the release notes.

Today I'll tackle imreconstruct, which performs morphological reconstruction. Reconstruction is a very useful operation that I've written about here before. For example, see my post from last year on opening by reconstruction. Several other Image Processing Toolbox functions call imreconstruct, including imclearborder, imfill, imhmax, imhmin, imextendedmax, imextendedmin, and imimposemin.

Let me use the opening by reconstruction example as a benchmark case.

url = 'http://blogs.mathworks.com/images/steve/2008/book_text.png';
text = imread(url);
imshow(text, 'InitialMagnification', 25)
title('918-by-2018 image displayed at 25% magnification')

The example task was to find letters containing vertical strokes by eroding with a vertical structuring element and then performing reconstruction.

se = strel(ones(51, 1));
marker = imerode(text, se);
text2 = imreconstruct(marker, text);
imshow(text2, 'InitialMagnification', 25)
title('Output image displayed at 25% magnification')

So how long does that call to imreconstruct take in R2009b? I'll use my function timeit, which you can download from the MATLAB Central File Exchange.

timeit(@() imreconstruct(marker, text))

ans =

    0.0241

That time is about 45 times faster than the same operation performed in R2009a. Note that I'm running on two-core computer; the improved imreconstruct is multithreaded, so the performance improvement would be greater on a four-core computer, for example.

Now let's time gray-scale reconstruction. I'll make a 1024-by-1024 test image and compute a marker image by subtraction. This kind of operation is often used to suppress small peaks in an image.

I = repmat(imread('rice.png'), 4, 4);
marker = I - 20;
timeit(@() imreconstruct(marker, I))

ans =

    0.0201

This time is about 30 times faster than R2009a, again running on my two-core laptop.

Now for some key caveats you should know. For now, the performance improvements described here only work for 2-D inputs that are uint8, uint16, or single, and only when the specified connectivity is 4 or 8. We'll be working in the future to extend the speed improvements to other inputs.

Get the MATLAB code

Published with MATLAB® 7.9