Steve on Image Processing

Concepts, algorithms & MATLAB

Aliasing and image resizing – part 3 4

Posted by Steve Eddins,

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.

  1. 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).
  2. 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.

  1. Convert the discrete-domain signal, $f[n]$, to a continuous-domain signal, $f(x)$.
  2. 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)$.
  3. 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

4 CommentsOldest to Newest

Cris replied on : 1 of 4

Steve,

Nice series on aliasing!

I think the key here is to understand that the two low-pass filters (interpolation filter of step 1 and blurring filter of step 2) have the same shape but different width. When multiplied in the Fourier domain, they lead to a single low-pass filter that is the narrower of the two (wider in the spatial domain). Thus, if filter #2 is wider, you use that one, otherwise you use #1. This is why you can simply make the interpolation filter wider when down-sampling, and don’t change anything when up-sampling.

Steve Eddins replied on : 2 of 4

Cris—Thanks, I appreciate the kind words! Yes, the way you formulated it is the way I learned it and the way I have taught it. I didn’t want to get too deep into signal processing theory, though. For example, I didn’t want to get into the difference between an ideal frequency response and the frequency responses of the interpolating kernels. That might be a topic for another day!

Jon replied on : 3 of 4

Nice series, thanks for explaining this in clear and simple way.
One question, if i scale 1/2 by averaging pixels and plot the interpolation kernel, would it just be a pulse at +-1 with amplitude of 1?
Thanks.

Steve Eddins replied on : 4 of 4

Jon—Yes. Here are a couple of details to keep in mind: (a) the pulse height would be 0.5, and (b) the pulse needs to be equal to 0.5 over a half-open interval. That is, -1 < = x < 1. Or -1 < x <= 1, depending on which way you want to break ties.

Add A Comment

What is 4 + 10?

Preview: hide