Steve on Image Processing

February 1st, 2010

Updates on old business: upslope area and unit testing

A couple of papers have been published recently that are related to some of my old blog topics.

Almost exactly a year ago I posted an xUnit-based unit test harness for MATLAB to the MATLAB Central File Exchange, and I blogged about it here.

Later in the year I wrote an article about it for Computing in Science and Engineering, a joint publication of the IEEE Computer Society and the American Institute of Physics. The article, "Automated Software Testing for MATLAB," appears in the November 2009 issue.

In 2007 I did a series of blog post on developing MATLAB algorithms for a topological measurement called upslope area, a measurement related to digital elevation models (DEMs). The upslope area of a point (pixel) is the total area on the terrain above that point that drains through it. I formulated the computation of upslope area as a sparse linear system, which let me put the MATLAB backslash operator to good use, and I put my upslope area functions on the File Exchange.

At about the same time, File Exchange contributor Wolfgang Schwanghart was developing a similar approach. His paper (with Nikolaus Kuhn), "TopoToolbox: A set of Matlab functions for topographic analysis," is in the January 2010 issue of Environmental Modelling & Software. TopoToolbox is available here.

January 18th, 2010

Relationship between continuous-time and discrete-time Fourier transforms

Previously in my Fourier transforms series I've talked about the continuous-time Fourier transform and the discrete-time Fourier transform. Today it's time to start talking about the relationship between these two.

Let's start with the idea of sampling a continuous-time signal, as shown in this graph:

Mathematically, the relationship between the discrete-time signal and the continuous-time signal is given by:

(When I write equations involving both continuous-time and discrete-time quantities, I will sometimes use a subscript "c" to distinguish them.)

The sampling frequency is (in Hz) or (in radians per second).

The discrete-time Fourier transform of is related to the continuous-time Fourier transform of as follows:

But what does that mean? There are two key pieces to this equation. The first is a scaling relationship between and : . This means that the sampling frequency in the continuous-time Fourier transform, , becomes the frequency in the discrete-time Fourier transform. The discrete-time frequency corresponds to half the sampling frequency, or .

The second key piece of the equation is that there are an infinite number of copies of spaced by .

Let's look at a graphical example. Suppose looks like this:

Note that equals zero for all frequencies . This is what we mean when we say a continuous-time signal is band-limited. The frequency is called the bandwidth of the signal.

The discrete-time Fourier transform of looks like this:

where . As I mentioned before, normally only one period of is shown:

For this example, then, between and looks just like a scaled version of .

Next time we'll consider what happens when doesn't look like . In other words, we're about to tackle aliasing.


Get the MATLAB code

Published with MATLAB® 7.9

January 11th, 2010

About the unused-argument syntax in R2009b

Last fall Loren wrote a blog post about a new syntax in R2009b for ignoring function inputs or function outputs. For example, suppose you call sort and only want the second output. In R2009b you can do that using this syntax:

 [~, idx] = sort(A);

The tilde character here is used as a placeholder to indicate that you don't need the first output argument of sort.

For another example, suppose you write a function that has to take three input arguments (because another function is always going to pass it three arguments), but your function doesn't need the second argument. Then you can write the first line of your function this way:

 function out = myfun(A, ~, C)

This new syntax has drawn a startling amount of discussion, both in the comments on Loren's post as well as in the comp.soft-sys.matlab newsgroup. (See the thread "Getting indexes of rows of matrix with more than n repetitions," for example.)

Responses to this new syntax have fallen into roughly four categories:

  • "Finally, I've been waiting for this!"
  • Clarification questions about how it works and how to use it.
  • Complaints about the specific syntax chosen and suggestions for alternatives.
  • "The intro of the tilde op was nothing but a big, useless blunder."

The passionate arguments in comp.soft-sys.matlab caught my eye and have prompted me to add my two cents here. Although I have no particular expectation of changing anyone's mind, I thought it might be interesting to address one particular question that was expressed well by Matt Fig: "If something already works, why complicate things?"

Matt means that we already have at least a couple of ways to ignore an output variable. The first is to use a variable name that (hopefully) makes the programmer's intent clear:

  [unused, idx] = sort(A);

The second technique takes advantage of the way MATLAB assigns function outputs from left to right:

  [idx, idx] = sort(A);

So why go to the trouble of introducing a new syntax?

Well, you'll probably get a somewhat different answer from every MATLAB developer you ask. Some might even agree. (I'll note, however, that this was one of the least-controversial MATLAB syntax proposals ever considered by the MATLAB language team.)

The proposed syntax was originally considered years ago, sometime around 2001. The proposal was approved internally at that time but then didn't get implemented right away because of competing priorities.

To understand why we eventually did decide to go ahead and implement it, you have to understand how this fairly minor syntactic issue is connected to our long-term efforts to do something much more important: provide automatic (and helpful!) advice to MATLAB users about their code. If you've used the MATLAB Editor at all over the last few years, you've probably observed how it marks code sections and makes suggestions to you. These suggestions fall into several categories, such as errors, possible errors, performance, new features, etc.

One of the things we try to flag is potential programming errors. A coding pattern that very often indicates a programming error is when you save a computed value by assigning it to a variable, but then you never use that saved value. At the very least, that pattern may indicate "dead code," or code that was doing something useful at one time but is now just cruft.

Both of the conventions mentioned above ([unused,idx] = sort(A) and [idx,idx] = sort(A)) exhibit this pattern of computing and saving a value and then never using it.

That might not be obvious for the [idx,idx] = sort(A) case, but it's true. The first output of sort is assigned to the variable idx. Then the second output is assigned to idx, causing the first output value to be discarded without ever being used.

These cases aren't programming errors, though, because we've given you no other way to ignore outputs. We can't automatically and reliably distinguish between the intentional [idx,idx] = sort(A) and the similar-looking [y,y] = foobar(x) that is the result of a typo.

I think it's important for the programmer to communicate his or her intent very clearly (which is why I tend to prefer the [unused,idx] = sort(A) convention). It is useful to have a way to communicate intent syntactically instead of by convention. And the new syntax helps us, in a small way, progress toward our long-term goal of helping MATLAB users write better MATLAB code.

OK, fire away!


Get the MATLAB code

Published with MATLAB® 7.9

January 8th, 2010

In my mailbox this week: MCALab, astronomy, dmperm, and Sudoku

Some interesting IEEE articles arrived in my mailbox this week. Computing in Science and Engineering has "MCALab: Reproducible Research in Signal and Image Decomposition and Inpainting," by Fadili, Starck, Elad, and Donoho. MCA means morphological component analysis, which is a new term for me. MCALab is an open-source library of MATLAB routines for signal and image decomposition and inpainting. The library is available from www.morphologicaldiversity.org/mcalab.html.

This month's Signal Processing Magazine has a special section on astronomy and cosmology packed with image processing theory and applications, including:

  • "Synthetic Aperture Radio Telescopes," by Levanda and Leshem
  • "Calibration Challenges for Future Radio Telescopes," by Wijnholds, van der Tol, Nijboer, and van der Veen
  • "Bayesian Source Separation for Cosmology," by Kuruoglu
  • "Precision Cosmology with the Cosmic Microwave Background," by Cardoso
  • "Cosmic Microwave Background Images," by Herranz and Vielva
  • "Light on Dark Matter with Weak Gravitational Lensing," by Pires, Starck, and Réfrégier.
  • "Multidimensional Image Reconstruction in Astronomy," by Kamalabadi
  • "Image Reconstruction in Optical Interferometry," by Thiébaut and Giovannelli
  • "Imaging with Linc-Nirvana," by Bertero, Boccacci, Desiderà, La Camera, Carbillet, and Lantéri

Finally, another article in Computing in Science and Engineering caught my eye: "Whip Until Solved," by Sullivan. This article describes a method for solving Sudoku puzzles that makes use of the MATLAB function dmperm. I confess that I can't think of a connection between Sudoku and image processing (take that as a challenge, dear reader), and I've never even solved a Sudoku puzzle. But I have written about dmperm, though, in my series on connected components labeling. And, coincidentally, MATLAB creator Cleve Moler wrote an article in the most recent MathWorks News & Notes about a different way to use MATLAB to solve Sudoku puzzles.

I've got an item on my to-do list that says "Study Cleve's Sudoku code to see how it works." Not sure when I'll get the time to do that, though.

December 31st, 2009

Discrete-time Fourier transform (DTFT)

In the last two posts in my Fourier transform series I discussed the continuous-time Fourier transform. Today I want to start getting "discrete" by introducing the discrete-time Fourier transform (DTFT).

The DTFT is defined by this pair of transform equations:

Here x[n] is a discrete sequence defined for all n:

I am following the notational convention (see Oppenheim and Schafer, Discrete-Time Signal Processing) of using brackets to distinguish between a discrete sequence and a continuous-time function. n is unitless. The frequency-domain variable, , is continuous with units of radians.

Note that is periodic with period .

Here are a few common transform pairs:

Unit Impulse

DTFT of Unit Impulse

Rectangular Pulse

DTFT of Rectangular Pulse

Note that the DTFT of a rectangular pulse is similar to but not exactly a sinc function. It resembles the sinc function between and , but recall that is periodic, unlike the sinc function.

Cosine

DTFT of Cosine

The DTFT of a discrete cosine function is a periodic train of impulses:

I updated the above plot on 6-Jan-2010 to show the location of the impulses. -SE

Because of the periodicity of it is very common when plotting the DTFT to plot it over the range of just one period: . For example, the DTFT of the rectangular pulse will most often be shown like this:

Next time I'll discuss the relationship between the continuous-time and the discrete-time Fourier transforms. Until then, Happy New Year everyone!


Get the MATLAB code

Published with MATLAB® 7.9

December 18th, 2009

Image Enhancement: A Medley

Can I post two silly-season items in one day? Hey, it's my blog, so why not!

One of my posts that received the most comments was the one about image processing in the movies.

Since so many readers commented on that post, I thought you all might enjoy this video medley of image enhancement in the movies and TV.

Remember, if your image processing program isn't working right, it's probably because "the eigenvalue is off."

A very different take on the Fourier transform

This is my silly season entry (from xkcd):

December 16th, 2009

Continuous-time Fourier transform of windowed cosine

Last week I showed a couple of continuous-time Fourier transform pairs (for a cosine and a rectangular pulse). Today I want to follow up by discussing one of the ways in which reality confounds our expectations and causes confusion. Specifically, when we're talking about real signals and systems, we never truly have an infinitely long signal.

To refresh your memory, here are the ideal cosine signal and its continuous-time Fourier transform plots again:

The dots at the left and right of the cosine plot are meant to remind you that the cosine signal is defined for all t.

But what if we modify the cosine signal so that we have only a finite number of periods? Let's start with a single period, for example.

% N = number of desired cosine periods
N = 1;

t1 = -N*pi;
t2 = N*pi;
t = linspace(t1, t2, 1000);
x = cos(t);
plot(t,x)
ylim([-1.25 1.25])
xlabel('t')
title('x(t) = cos(t)')

In signal processing jargon you'll see this kind of signal referred to as a windowed cosine with a rectangular window. That is, the signal is an infinite-extent cosine multiplied by a finite-extent rectangular pulse.

If I've done my math correctly, the continuous-time Fourier transform of a cosine signal truncated to a finite number of periods is:

where N is the number of periods.

(And if I haven't done my math correctly, I'm sure a sharp-eyed reader will let me know. I feel like I'm back at Georgia Tech working on homework sets. Except that I'm doing it during the day instead of late at night.)

omega = linspace(-5.5, 5.5, 500);
X = sin(pi*N*(omega-1))./(omega-1) + sin(pi*N*(omega+1))./(omega+1);
plot(omega, X)
title('X(\Omega)')
xlabel('\Omega (rad/s)')

Well, that doesn't look very much like the continuous-time Fourier transform for the ideal cosine! If you run this code yourself, try zooming in to look closely at those peaks. You'll see they aren't even located precisely at 1.0 rad/s.

Let's try 10 periods.

N = 10;
t1 = -N*pi;
t2 = N*pi;
t = linspace(t1, t2, 1000);
x = cos(t);
plot(t,x)
ylim([-1.25 1.25])
xlabel('t')
title('x(t) = cos(t)')
omega = linspace(-5.5, 5.5, 500);
X = sin(pi*N*(omega-1))./(omega-1) + sin(pi*N*(omega+1))./(omega+1);
plot(omega, X)
title('X(\Omega)')
xlabel('\Omega (rad/s)')

That's closer, but there's still a lot of funny stuff going on. The more periods we use, the closer the Fourier transform gets to the ideal case, but for any finite number of periods we're going to get this decaying ripple effect.

Later, when we start constructing synthetic spatial sinusoid images (like in this post), and what we see in the frequency domain looks "messy," this is one of the big reasons. All our images have finite extent! We'll be coming back to this idea in future posts.

Have you heard the term bandlimited? It refers to a signal that has no frequency energy above a certain cut-off frequency. Ponder this notion: No finite-extent signal is bandlimited. That is, any signal that is nonzero over only a finite portion of the time domain has nonzero frequency energy across the entire frequency domain.

I think that next time I'll be ready to start talking about the discrete-time Fourier transform, or DTFT. Conceptually we are traveling methodically toward the discrete Fourier transform, or DFT, which is what the MATLAB function fft computes. Along the way we'll figure out how all three forms (continuous-time Fourier transform, discrete-time Fourier transform, and discrete Fourier transform) relate to each other.


Get the MATLAB code

Published with MATLAB® 7.9

December 11th, 2009

Continuous-time Fourier transform basics

In my planned (well, partially planned) discussion on Fourier transforms, I'll focus on three of the four types I listed in my November 23 post:

  • Continuous-time Fourier transform
  • Discrete-time Fourier transform
  • Discrete Fourier transform

The existence of multiple transform flavors, as well as the details of their relationships, is at the heart of much of the confusion on this topic.

Let's start with the continuous-time Fourier transform. (When the context makes it clear whether I'm talking about the continuous-time or the discrete-time flavor, I'll often just use the term Fourier transform.)

The continuous-time Fourier transform is defined by this pair of equations:

There are various issues of convention and notation in these equations:

  • You may see a different letter used for the frequency domain ( or f, for example). I am in the habit of using for the continuous-time Fourier transform and for the discrete-time Fourier transform.
  • You may see i instead of j used to represent . I tend to follow the electrical engineering tradition of using j.
  • You may see terms appearing in the exponent of e and not in front of the inverse transform integral.
  • You may see the signs of the exponent terms switched in the transform equations (that is, the minus sign in the exponent appears in the inverse transform instead of the forward transform).

All these variations exist because we don't already have enough to worry about in the rest of our lives.

With the equations I use, the frequency domain unit is angular frequency (radians/second).

There are so many useful Fourier transform properties and transform pairs that it's hard for to pick the bare minimum necessary for the ideas I want to convey. Today I'll just show you two of the most essential Fourier transform pairs in signal processing applications.

Here's a cosine signal:

Plots corrected December 14 thanks to help from Mark Andrews.

And here is its Fourier transform:

This is what most people who have some knowledge of the Fourier transform expect to see. A signal containing a single frequency (here the frequency is 1 rad/s) has all its frequency domain energy concentrated at that single frequency.

The second pair is a rectangular pulse in the time domain and a sinc function in the frequency domain.

I'd like to call your attention especially to the dots ("...") at the left and right on the cosine and sinc function plots. The dots are there to remind you that these functions have infinite extent. That is, they are nonzero over the entire domain.

That's an important thing to keep in mind. I'll come back to that point next time.

One final note: I've started using the category "Fourier transforms" for posts on this topic. You can see all the related posts by clicking on the category link on the right side of the page.


Get the MATLAB code

Published with MATLAB® 7.9

December 7th, 2009

Example of Fourier terminology confusion

This exchange of comments between Cris and me (see comments 11-16) is a good example of the kind of terminology confusion that surrounds Fourier transforms. Cris and I both have a background in signal processing theory, but we were using the term "DFT" to refer to two different transform types, and it took a little while for us to figure that out.

Steve Eddins manages the Image & Geospatial development team at The MathWorks and coauthored Digital Image Processing Using MATLAB. He writes here about image processing concepts, algorithm implementations, and MATLAB.

  • Steve: Dansheng—The following code will give you the binary image corresponding to a single object. L =...
  • Dansheng Song: Hi Steve, Your posts are very helpful for Matlab users. Do you know how can I get the binary image by...
  • DP: Thanks Steve. it did work. now can easily calculate DSC.
  • Steve: DP— A & B
  • DP: Hi Steve, Dice similarity coefficient is defined to find the overlapped regions between two objects in an image....
  • Steve: Searching— axis on
  • Steve: DP—Sorry, but I’ve never heard of dice similarity coefficients.
  • searching: Hey, I found this place extremely useful. I have a question though. Although I manage to plot my points on...
  • DP: Hi Steve, thanks for your email. i could not seperate those connected objects. Instead i am trying to compare...
  • Gopesh: Excellent code..was very helpful to me to change the background color of mri images and thereby improving the...

These postings are the author's and don't necessarily represent the opinions of The MathWorks.