With the R2015a release a couple of years ago, the Image Processing Toolbox added the function imgaussfilt. This function performs 2-D Gaussian filtering on images. It features a heuristic that automatically switches between a spatial-domain implementation and a frequency-domain implementation.... read more >>

]]>With the R2015a release a couple of years ago, the Image Processing Toolbox added the function `imgaussfilt`. This function performs 2-D Gaussian filtering on images. It features a heuristic that automatically switches between a spatial-domain implementation and a frequency-domain implementation.

I wanted to check out the heuristic and see how well it works on my own computer (a 2015 MacBook Pro).

You can use `imgaussfilt` this way:

```
rgb = imread('peppers.png');
sigma = 10;
rgb_smoothed = imgaussfilt(rgb,sigma);
imshow(rgb)
```

```
imshow(rgb_smoothed)
title('Gaussian smoothed, \sigma = 10')
```

With this call, `imgaussfilt` automatically chooses between a spatial-domain or a frequency-domain implementation. But you can also tell `imgaussfilt` which implementation to use.

rgb_smoothed_s = imgaussfilt(rgb,sigma,'FilterDomain','spatial'); rgb_smoothed_f = imgaussfilt(rgb,sigma,'FilterDomain','frequency');

I'll use this syntax, together with the `timeit` function, to get rough timing curves for the two implementation methods. The heuristic used by `imgaussfilt` uses a few different factors to decide, including image size, Gaussian kernel size, single or double precision, and the availability of processor-specific optimizations. For my computer, with a 2000-by-2000 image array, the cross-over point is at about $\sigma = 50$. That is, `imgaussfilt` uses the spatial-domain method for $\sigma < 50$, and it uses a frequency-domain method for $\sigma > 50$.

Let's check that out with a set of $\sigma$ values ranging from 5 to 200.

I = rand(2000,2000); ww = 5:5:200; for k = 1:length(ww) t_s(k) = timeit(@() imgaussfilt(I,ww(k),'FilterDomain','spatial')); end for k = 1:length(ww) t_f(k) = timeit(@() imgaussfilt(I,ww(k),'FilterDomain','frequency')); end plot(ww,t_s) hold on plot(ww,t_f) hold off legend('spatial','frequency') xlabel('\sigma') ylabel('time (s)') title('imgaussfilt - 2000x2000 image')

As you can see, the cross-over point of the timing curves is about the same as the threshold used by `imgaussfilt`. As computing hardware and optimized computation libraries evolve over time, though, the ideal cross-over point might change. The Image Processing Toolbox team will need to check this every few years and adjust the heuristic as needed.

Get
the MATLAB code

Published with MATLAB® R2016b

Today I'll try to wrap up my discussion about how aliasing affects image resizing and about how the imresize function tries to prevent it. (This is called antialiasing.) Let me show you where we are going. Here is the zone plate image I showed last time.... read more >>

]]>Today I'll try to wrap up my discussion about how aliasing affects image resizing and about how the `imresize` function tries to prevent it. (This is called *antialiasing*.) Let me show you where we are going. Here is the zone plate image I showed last time.

```
Z = imzoneplate(501);
imshow(Z)
title('Z')
```

Here's what happens when we shrink `Z` by throwing away samples.

```
Z4 = Z(1:4:end,1:4:end);
imshow(Z4)
title('Z4')
```

And here's what `imresize` does (with the default behavior):

```
Z4_imresize = imresize(Z,0.25);
imshow(Z4_imresize)
title('imresize(Z,0.25)')
```

Before I get started describing the `imresize` approach to antialiasing, I should say that there are several algorithm approaches to antialiasing in image resizing. This is just one. It is based on classical ideas from digital signal processing.

I'll use the notation that $f[n]$ is a function of a discrete variable $n$, and $f(x)$ is a function of a continuous variable $x$. Here's a way to think about the image resizing problem that is theoretically useful. To simplify the discussion for a bit, I'll talk in terms of 1-D signals. (Warning: people with a graduate-level background in digital signal processing will probably find this explanation too hand-wavy. Everybody else will find it too jargon-y. I'm doing my best.)

Suppose we want to convert a 1-D signal, $f[n]$, to another 1-D signal that has a different sampling rate.

- Convert the discrete-domain signal, $f[n]$, to a continuous-domain signal, $f(x)$. Theoretically, this step involves the use of a continuous-domain filter, $h(x)$. You can think of this step as a form of
**interpolation**(for some choices of $h(x)$, anyway). - Convert the continuous-domain signal, $f(x)$, to a new discrete-domain signal, $g[n]$, by sampling $f(x)$ at the desired rate: $g[n] = f(nT)$.

There's a problem with the two-step procedure above. If the sampling rate for $g[n]$ is reduced from the original sampling rate of $f[n]$, then the procedure is susceptible to aliasing distortion if $f[n]$ contains high frequencies that cannot be represented at the lower sampling rate. To address this problem, let's insert another step in the procedure.

- Convert the discrete-domain signal, $f[n]$, to a continuous-domain signal, $f(x)$.
- Pass $f(x)$ through a lowpass filter to remove high-frequency components that would cause aliasing distortion at the desired sampling rate. Let's call this filter $h_a(x)$.
- Convert the continuous-domain signal, $f(x)$, to a new discrete-domain signal, $g[n]$, by sampling $f(x)$ at the desired rate: $g[n] = f(nT)$.

One key to understanding the `imresize` approach is to recognize that filtering with two filters in succession, in this case $h(x)$ and $h_a(x)$, can be replaced by one filter whose frequency response is the product of the original two filters. The second key to understanding the `imresize` algorithm is that the desired frequency response for the antialiasing filter, $h_a(x)$, depends how much we are changing the sampling rate. If we shrink the original signal more, then we have to remove more of the high frequencies with $h_a(x)$.

To state that more compactly (while waving my hands fairly vigorously), you can use just **one** continuous-domain filter to simultaneously accomplish interpolation and antialiasing.

Let me show you how this works for the cubic interpolation kernel that `imresize` uses by default. Here is the code for the cubic interpolation kernel:

function f = cubic(x) % See Keys, "Cubic Convolution Interpolation for Digital Image % Processing," IEEE Transactions on Acoustics, Speech, and Signal % Processing, Vol. ASSP-29, No. 6, December 1981, p. 1155.

absx = abs(x); absx2 = absx.^2; absx3 = absx.^3;

f = (1.5*absx3 - 2.5*absx2 + 1) .* (absx <= 1) + ... (-0.5*absx3 + 2.5*absx2 - 4*absx + 2) .* ... ((1 < absx) & (absx <= 2));

(This code is inside the file imresize.m.) And here's what the interpolation kernel looks like:

```
figure
fplot(@cubic,[-2.5 2.5],'LineWidth',2.0)
```

The function `imresize` then modifies the interpolation by stretching it out. The stretch factor depends directly how much we are shrinking the original signal. Here's the `imresize` code that modifies the interpolation kernel:

if (scale < 1) && (antialiasing) % Use a modified kernel to simultaneously interpolate and % antialias. h = @(x) scale * kernel(scale * x); kernel_width = kernel_width / scale; else % No antialiasing; use unmodified kernel. h = kernel; end

(I don't always comment my code well, but I think I did OK here.)

Notice that we don't modify the kernel at all when we are **growing** a signal (`scale > 1`) instead of shrinking it (`scale < 1`).

Here are three versions of the cubic interpolation kernel: the original one, one modified for shrinking by a factor of 2, and one modified by shrinking by a factor of 4.

h1 = @cubic; h2 = @(x) 0.5 * h1(0.5*x); h4 = @(x) 0.25 * h1(0.25*x); fplot(h1,[-10 10],'LineWidth',2) hold on fplot(h2,[-10 10],'LineWidth',2) fplot(h4,[-10 10],'LineWidth',2) hold off legend('h_1(x)','h_2(x)','h_4(x)')

When we spread out the interpolation kernel in this way, the interpolation computation averages pixel values from a wider neighborhood to produce each output pixel. That is additional smoothing, which gives us the desired antialiasing effect.

Here is the `imresize` output for the zone plate image again.

imshow(Z4_imresize)

Only the rings near the center are visible. That's because those rings had lower spatial frequency in the original, and they can be successfully represented using only one-fourth the sampling rate. The higher-frequency rings further away from the center have been smoothed away into a solid gray.

You can experiment yourself with `imresize` by disabling antialiasing. Here's how to shrink an image with antialiasing turned off:

```
Z4_imresize_noaa = imresize(Z,0.25,'Antialiasing',false);
imshow(Z4_imresize_noaa)
```

To give credit where credit is due, the algorithm used by `imresize` was inspired by the article "General Filtered Image Rescaling," by Dale Schumacher, in *Graphics Gems III*, Morgan Kaufmann, 1994.

Do you have other questions about `imresize`? Let me know by leaving a comment.

Get
the MATLAB code

Published with MATLAB® R2016b

In my 03-Jan-2017 post, I introduced my new topic: aliasing, especially as it relates to resizing images. Today I will get more specific about what aliasing actually is.... read more >>

]]>In my 03-Jan-2017 post, I introduced my new topic: *aliasing*, especially as it relates to resizing images. Today I will get more specific about what aliasing actually is.

I'll use a 1-D sequence to illustrate. Here is part of a 1-D sequence with a period of 8.

a = repmat([1 1 1 1 0 0 0 0],1,9)

a = Columns 1 through 13 1 1 1 1 0 0 0 0 1 1 1 1 0 Columns 14 through 26 0 0 0 1 1 1 1 0 0 0 0 1 1 Columns 27 through 39 1 1 0 0 0 0 1 1 1 1 0 0 0 Columns 40 through 52 0 1 1 1 1 0 0 0 0 1 1 1 1 Columns 53 through 65 0 0 0 0 1 1 1 1 0 0 0 0 1 Columns 66 through 72 1 1 1 0 0 0 0

Now let's "shrink" the sequence by sampling it with a spacing that is greater than 1, using nearest-neighbor interpolation. First, let's try a spacing of 4/3.

f = 4/3; a_1 = a(round(1:f:end))

a_1 = Columns 1 through 13 1 1 1 0 0 0 1 1 1 0 0 0 1 Columns 14 through 26 1 1 0 0 0 1 1 1 0 0 0 1 1 Columns 27 through 39 1 0 0 0 1 1 1 0 0 0 1 1 1 Columns 40 through 52 0 0 0 1 1 1 0 0 0 1 1 1 0 Columns 53 through 54 0 0

The shrunken sequence has a period of 6. That's consistent with the original sequence period and the shrink factor: $8/(4/3) = 6$.

Now let's increase the shrink factor to 1.6.

f = 1.6; a_2 = a(round(1:f:end))

a_2 = Columns 1 through 13 1 1 1 0 0 1 1 1 0 0 1 1 1 Columns 14 through 26 0 0 1 1 1 0 0 1 1 1 0 0 1 Columns 27 through 39 1 1 0 0 1 1 1 0 0 1 1 1 0 Columns 40 through 45 0 1 1 1 0 0

The sequence `a_2` has a period of 5. Again, that is consistent with the original sequence period and shrink factor: $8/1.6 = 5$.

Let's bump the shrink factor up to 4.

f = 4; a_3 = a(round(1:f:end))

a_3 = Columns 1 through 13 1 0 1 0 1 0 1 0 1 0 1 0 1 Columns 14 through 18 0 1 0 1 0

The sequence `a_4` has a period of 2, which is still understandably straightforward: $8/4 = 2$. Note that 2 is the shortest period possible for a nonconstant periodic discrete sequence.

It gets more fun (and less straightforward) when the shrink factor goes above 4.

f = 6; a_4 = a(round(1:f:end))

a_4 = 1 0 0 1 1 0 0 1 1 0 0 1

Now the shrunken sequence's periodicity appears to be 5, and that's harder to understand because $8/6 \neq 5$.

To recap what we've seen so far: the period of the shrunken sequence decreased steadily as we increased the shrink factor, but only up to a point. Once the shrink factor increased above 4 (half the original period), the period of the shrunken sequence stopped decreasing and starting changing in a different fashion.

This is an example of what signal processing folk call *aliasing*. The term is used in the sense of one frequency masquerading as another.

This happens when the shrink factor becomes too high to preserve the period of the original sequence.

Now let me illustrate aliasing using a zone plate image. (You can get imzoneplate from the MATLAB Central File Exchange.)

```
Z = imzoneplate(501);
imshow(Z)
title('Z')
```

```
Z4 = Z(1:4:end,1:4:end);
imshow(Z4)
title('Z4')
```

That result hardly looks anything like the original image. The visual artifacts are caused by aliasing. The original image has high spatial frequencies near the boundaries that can't be represented when we shrink the image. Those frequencies, instead of just disappearing, get "aliased" into a different set of spatial frequencies.

In general, aliasing distortion can't be removed once it has happened. Instead, antialiasing techniques must be applied as part of the image resizing process. Next time, I'll discuss the particular antialiasing procedure used by `imresize`.

Get
the MATLAB code

Published with MATLAB® R2016b

Happy New Year!... read more >>

]]>Happy New Year!

I'm going to kick off 2017 with two or three posts about *antialiasing*. Specifically, I will discuss what antialiasing means in the context of image resizing, and I will explain how the `imresize` function does it.

Let me show you two pairs of images. Image A is a "zone plate" image that I made using my function `imzoneplate`. (You can get this function from the MATLAB Central File Exchange).

Image A1 shows a bunch of concentrating cirles, some with a bright and some with a dark center, arranged on a grid.

**Image A**

**Image A1**

OK, here are images B and B1. Image B is a texture image drawn from a photograph of blue jeans.

Image B1 shows a thumbnail image with another texture that is rotated roughly 90 degrees from the texture in image B.

**Image B**

**Image B1**

As you may have guessed, image A1 is a thumbnail-resized version of image A, and image B1 is a thumbnail-resized version of image B.

Here's a GIF animation that I made. It shows several different thumbnails of image B that are made with slightly different sizes. As the thumbnail sizes change, you can see the geometrical pattern actually seem to swing around in different directions.

Next time I'll talk about what causes this effect and what can be done about it.

Get
the MATLAB code

Published with MATLAB® R2016b

Earlier this summer, I was writing some color-space conversion code. At one point in the code, I had a Px3 matrix called `RGB`, which contained P colors, one per row. I also had a 1x3 vector, `v`. I needed to multiply each column of `RGB` by the corresponding element of `v`, like this:

RGB_c = [RGB(:,1)*v(1) RGB(:,2)*v(2) RGB(:,3)*v(3)];

But since I was using an internal developer build of MATLAB R2016 (released on September 14), I didn't type the code above. Instead, I typed this:

RGB_c = RGB .* v;

In R2016a and older MATLAB releases, that line of code produced an error:

>> RGB_c = RGB .* v Error using .* Matrix dimensions must agree.

In the new release, though, MATLAB **implicitly expands** the vector `v` to be the same size as the matrix `RGB` and then carries out the elementwise multiplication. I say "implicitly" because MATLAB does not actually make an in-memory copy of the expanded vector.

To read all about this change in MATLAB matrix arithmetic, head over to Loren's Art of MATLAB blog and read my guest post about it.

Get
the MATLAB code

Published with MATLAB® R2016b

Intrepid MathWorks application engineer Brett Shoelson recently got a user question that caught my attention. Consider an image containing text characters in outline form, such as this:... read more >>

]]>Intrepid MathWorks application engineer Brett Shoelson recently got a user question that caught my attention. Consider an image containing text characters in outline form, such as this:

```
url = 'http://blogs.mathworks.com/steve/files/MathWorks-address-binary.png';
bw = imread(url);
imshow(bw)
```

How can we fill in the text characters from their outlines without filling in the internal holes? If we just use `imfill` with the `'holes'` option, you can see that it doesn't give us the desired result.

bw_filled = imfill(bw,'holes'); imshow(bw_filled) title('Original with holes filled')

When I saw this problem, I thought that some combination of `imfill`, `imclearborder`, and logical operators could possibly solve it.

You've already seen `imfill`. Here's how `imclearborder` works.

url_sample = 'http://blogs.mathworks.com/images/steve/168/aug31.png'; bw_sample = imread(url_sample); imshow(bw_sample) title('imclearborder demonstration - input image')

```
bw_sample_clearborder = imclearborder(bw_sample);
imshow(bw_sample_clearborder)
title('imclearborder demonstration - output image')
```

You can see that any connected component touching any image border has been removed.

Going back to our task, I'm going to proceed first by identifying the background pixels that are inside the text characters. Roughly speaking, I'll tackle this phase by working "from the outside in."

First, let's identify and remove the pixels that are external to the characters.

```
bw2 = ~bw;
imshow(bw2)
title('Complement of original image')
```

```
bw3 = imclearborder(bw2);
imshow(bw3)
title('External pixels removed from the foreground')
```

Now let's complement the image and clear the borders again.

```
bw4 = ~bw3;
imshow(bw4)
title('Complement of bw3')
```

```
bw5 = imclearborder(bw4);
imshow(bw5)
title('After second border-clearing')
```

If we fill the holes in the image `bw5` above, and then take the exclusive-or of the result with `bw5` we'll be left only with the internal hole pixels inside the characters.

bw6 = imfill(bw5,'holes'); bw7 = xor(bw6,bw5); imshow(bw7) title('Internal hole pixels')

We're almost there. We can now use `bw6` to "fix up" the initial filled result, `bw_filled`, using an exclusive-or operation.

```
bw_final = xor(bw_filled,bw7);
imshow(bw_final)
title('Presto!')
```

Readers, how would you solve this problem? Brett and I think there might be a couple of other reasonable approaches. Let us know in the comments.

Get
the MATLAB code

Published with MATLAB® R2016b

Because I work regularly with developers on the MATLAB Math team, I see the questions that come to them from tech support. Today there was an interesting one.

A user observed that computing the `sin` function takes longer for very big numbers. The user asked, "Why is that?"

Well, first let's see if we can reproduce the observation. Start by making a vector with one million numbers in the interval [0,1].

rng(42) x = rand(1000000,1);

(I seeded the random number generator using `rng` so that the numbers below don't change each time I run this script.)

How long does it take to compute the sin of all those numbers?

f = @() sin(x); timeit(f)

ans = 0.002102730998500

Now repeat the timing experiment with another vector containing one million numbers in the interval [100000000,100000000+1].

x2 = x + 100000000; g = @() sin(x2); timeit(g)

ans = 0.027343141998500

The user is right, it does take longer!

Here's the straightforward explanation: For large numbers, it takes more computation to produce an accurate result.

To illustrate, let's look at the output of `sin` for one large number in MATLAB and compare it the result from another application. I'm going to use Excel for this comparison, but the results will be similar for code you write in C that calls the standard C math library function, `sin()`.

z = x2(1)

z = 1.000000003745401e+08

matlab_result = sin(z)

matlab_result = 0.734111589042897

excel_result = 0.734111669990523

excel_result = 0.734111669990523

These two results agree to only 6 digits!

Let's use the Symbolic Math Toolbox to compute `sin(z)` with 20 decimal digits of accuracy.

sym_result = sin(sym(z))

sym_result = sin(3355443212567481/33554432)

vpa(sym_result,20)

ans = 0.73411158904289725019

From this calculation, you can see that the MATLAB result is accurate to 15 decimal digits, whereas the Excel result is accurate to only 6 digits.

That's called sweating the details.

Get
the MATLAB code

Published with MATLAB® R2016b

This year also marks my tenth year of blogging. I still remember in 2005 when my friend Ned Gulley asked me, “Hey Steve, what do you think about writing a blog for MATLAB Central?” Two thoughts passed immediately through my head.

- That sounds like a really interesting opportunity, and it would be a nice way to push beyond my comfort zone.
- If you do your usual thing of analyzing and pondering Ned’s question for a few days, you’re going to get cold feet.

So I just said, “Sure, I’ll do it!” A few months later, in January 2006, a little bit after Loren started her Art of MATLAB blog, I made my first post.

That first year, I posted about algorithms (such as polygon scan conversion), signal processing concepts (such as separable convolution), Image Processing Toolbox tips and tricks (such as labeling labeled objects), and MATLAB in general (such as my series on how MATLAB image display works).

And that’s pretty much the way it’s worked ever since.

I recall that it took me almost four years to get up the nerve to start posting about Fourier transform concepts. It is such a misunderstood area, and people with different technical backgrounds learn the concepts in different orders and with different terms. But I finally dipped my toes in those waters in November 2009, and now there are 21 posts in the Fourier transforms category.

I’ll leave you with pointers to two other memorable posts. First, just before a summer holiday weekend in 2007, I invited readers to post their favorite uses and abuses of image processing in the movies. I loved the comment thread that resulted.

The other was in October 2006, when I broke my nine-year silence to explain the story behind the default MATLAB image.

That’s it for this post - my 500th!

]]>Ned commented that MATLAB Central was “giving a modern web-based home to a community that was already thriving.” That reminded me of how important MathWorks participation in that community has been to me personally.

Many, many moons ago, as the calendar changed and 1992 became 1993, I was a fairly unhappy electrical engineering professor. I was thinking about changing my career path, but I didn’t have any solid ideas about what to do outside academia.

I was a heavy MATLAB user, so I immediately subscribed when the new Usenet group, comp.soft-sys.matlab, launched in January. I soon noticed when Cleve Moler, creator of the original MATLAB, posted a note from MathWorks:

Hi. This is the first posting from The MathWorks to comp.soft-sys.matlab. I’m posting it from Stanford because we’re having trouble posting to newly establish groups, including comp.soft-sys.matlab, from the MathWorks machines in Massachusetts. We can post to old, established groups, but not to some new ones, including what is now the most important group for us. If anybody else is having similar difficulties, please let us know. Thanks to Michael Maurer at Stanford for initiating the group. The size of the vote in favor, and the level of participation this first week, are certainly good signs for an active, useful group.

At the MathWorks, we intend to follow the group’s activities, to respond individually to technical and informational queries, and to post general responses and announcements when we believe they are of broad interest. We certainly intend to respect the Net’s noncommercial culture.

We hope that many of the queries will be answered by other group participants. The discussion we’ve seen so far of MATLAB/Unix memory management is a good example of how this ought to work.

Thanks to all of you for your support of MATLAB. With your help, comp.soft-sys.matlab can be a valuable, and enjoyable, activity for all of us.

– Cleve Moler

Chairman and Chief Scientist

The MathWorks

And Cleve kept posting, always eager to teach people how to get their mathemetical computations done accurately and efficiently. For example, this 30-Jan-1993 post was in response to a question about computing factorials in MATLAB:

MATLAB does not have a factorial function, so n! must be computed by

`prod(1:n)`

or, possibly with roundoff error, by

`gamma(n+1)`

But, be careful, either of these quantities overflow for n > 170. If you actually want n! to compute something like

`n!/(k!*(n-k)!)`

for large n, you should use

`exp(sum(log(1:n))-sum(log(1:k))-sum(log(1:(n-k))))`

It’s not as slow as it looks, and it doesn’t overflow.

– Cleve Moler

Then there was this fellow named John Little, who started posting on the same day as Cleve and who kept up a steady stream of technical MATLAB nuggets like this one:

The MEX-file facility in MATLAB 4.0 includes a new routine called mexPutMatrix that allows you to plop a matrix directly back into the MATLAB workspace.

John Little

(At MathWorks, we generally call him “Jack.” He is our CEO.)

I noticed other MathWorks people jumping in, too. Clay Thompson answered one of my own questions in early February 1993. And then there was Loren Shure, who was answering a question about the `quad` function in March 1993. (Many of you know Loren from her Art of MATLAB blog.)

By following these MathWorks interactions in the early MATLAB community, I began to form a very favorable impression of the company and the people who worked there. When I saw Loren at a signal processing conference that spring, I introduced myself, and I volunteered to be a beta tester for the new image processing product that Loren could only hint about. By the fall of 1993, I had participated in the beta test and then interviewed (successfully) for a job. My family and I moved from Chicago to Massachusetts in December, and we’ve been here ever since!

]]>