On April 20, I wrote about an algorithm for solving the Eight Queens Problem. In this problem, the task is to place eight queens on a chessboard so that none of the queens is attacking any other queen. Here is one solution.... read more >>

]]>On April 20, I wrote about an algorithm for solving the Eight Queens Problem. In this problem, the task is to place eight queens on a chessboard so that none of the queens is attacking any other queen. Here is one solution.

After I published that post, I became curious to see what others might have submitted about this problem on the File Exchange. A little searching brought me to this submission by uu tii. When I looked at the submission, I was surprised to how little code it contained. This is it:

function M = eightqueens() N=8; T=N-2; rows=1:N; cols=perms(rows); S=size(cols,1); M=zeros(N,N,S); linearInd = sub2ind(size(M), repmat(rows',1,S), cols', repmat(1:S,N,1)); M(linearInd) = 1; dv=arrayfun(@(k)max([arrayfun(@(x)sum(diag(M(:,:,k),x)),-T:T),arrayfun(@(x)sum(diag(rot90(M(:,:,k)),x)),-T:T)]),1:S); M(:,:,dv>1)=[];

When I ran it, it returned an 8x8x92 array. The size caught my eye because there are exactly 92 solutions to the Eight Queens Problem (including rotations and flips).

```
M = eightqueens;
whos M
```

Name Size Bytes Class Attributes M 8x8x92 47104 double

Let's peek at one of the planes of `M`.

M(:,:,45)

ans = 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0

Indeed, that is a solution. So, what's going on in this small chunk of code? Well, as I see it, this code has three conceptual chunks. The first chunk, consisting of eight lines, is the longest. Let's run it.

```
N=8;
T=N-2;
rows=1:N;
cols=perms(rows);
S=size(cols,1);
M=zeros(N,N,S);
linearInd = sub2ind(size(M), repmat(rows',1,S), cols', repmat(1:S,N,1));
M(linearInd) = 1;
whos M
```

Name Size Bytes Class Attributes M 8x8x40320 20643840 double

This seems to be a collection of 40320 boards. There aren't that many solutions, so the collection must contain many boards that aren't valid solutions. For example:

M(:,:,45)

ans = 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0

There are four queens lined up on the main antidiagonal, so this board is not a solution.

This set of boards has been constructed by making the observation that any solution must have exactly one queen on each row and exactly one queen on each column. The array `M` contains every board that has this property. All the computational magic here is in the line `cols=perms(rows)`. The rest of this chunk is some indexing computations to initialize and fill in the array `M` with queens at the locations given by `rows` and `cols`.

The second chunk is this line:

dv=arrayfun(@(k)max([arrayfun(@(x)sum(diag(M(:,:,k),x)),-T:T),arrayfun(@(x)sum(diag(rot90(M(:,:,k)),x)),-T:T)]),1:S);

Wow, there's a lot going on there. I would summarize it this way: For every board `k` in M, sum every diagonal and every antidiagonal and take the maximum of those sums. The idea is that `M` was constructed so that every board has exactly one queen on each row and exactly one queen on each column, so you really only have to check the diagonals of a board to determine whether it is a solution. Let's look at `dv` for the 45th board, which I displayed above.

dv(45)

ans = 4

Remember my comment above that board 45 has 4 queens lined up on the main antidiagonal? That's why `dv(45)` is 4.

The last conceptual chunk is eliminate all the boards in the `M` that are not solutions. We do that by eliminating all boards for which `dv > 1`.

```
M(:,:,dv>1)=[];
whos M
```

Name Size Bytes Class Attributes M 8x8x92 47104 double

And there it is. That's a very clever use of `perms` and indexing and `arrayfun` to generate all possible solutions. Thanks, uu tii.

Get
the MATLAB code

Published with MATLAB® R2017a

Today I want to tackle a classic algorithm puzzle known as the Eight Queens Problem. In this problem, your task is to arrange eight queens on a chessboard so that none of the queens is attacking any of the others. Here is one solution.... read more >>

]]>Today I want to tackle a classic algorithm puzzle known as the *Eight Queens Problem*. In this problem, your task is to arrange eight queens on a chessboard so that none of the queens is attacking any of the others. Here is one solution.

The problem is often restated as the *N Queens Problem*, which is placing N queens on an N-by-N board.

This is more of a general algorithm topic than an image processing topic, although some of the concepts and MATLAB techniques that I show might be useful to people doing image processing work. For some reason, this problem has been on my mind recently, so I decided to play with it a bit in MATLAB. I ended up with both a working algorithm and an interesting visualization. Here is the result for a 6x6 chessboard.

Solving the eight queens problem was the first "real" and nontrivial program that I can recall writing. I wrote it in Basic for a math class project sometime around 1981. My school didn't actually have a computer at that point; I used a teletype terminal (no monitor!) and a dial-up modem to communicate with the county school system's mainframe, which was located some miles away. As I recall, the program used a pure brute-force solution, generating boards randomly and testing each one to see if it was a valid solution. I knew almost nothing about algorithms then.

When I revisited the problem this week, I used a *recursive divide-and-conquer technique with backtracking*. I'll explain that in two parts, starting with the *recursive divide-and-conquer* part. The most famous algorithm in signal processing, the fast Fourier transform (FFT) is based on this idea. The discrete Fourier transform of an $N$-point sequence can be written in terms of the discrete Fourier transforms of two different subsets of the original sequence.

Here is one way to use divide-and-conquer on the N queens problem. We'll write a routine that places N queens on an N-column board in two steps:

A. Place one queen on a safe square somewhere on the first column of the board.

B. Call the routine recursively to place N-1 queens on an (N-1)-column board.

The *backtracking* part is necessary because, unlike when computing the FFT, one of the substeps A or B might fail.

If step B fails, then we *backtrack*. That is, we undo step A, find a different safe square on the first column of the board, and then do the recursive step B call again. At any time, if there are no remaining safe squares on the first column to try, then the solution attempt has failed and the routine returns.

Here's some MATLAB code to implement this idea.

function [board,succeeded] = placeQueens(board,col)

N = size(board,1); if col > N % There are no columns left. All the queens % have been placed. Return with success. succeeded = true; return end

safe = false(1,N); for row = 1:N safe(row) = isSafe(board,row,col); end

if ~any(safe) % There are no safe squares on this column. % Return with failure. succeeded = false; return end

for row = 1:N if safe(row) board(row,col) = 1; [board,succeeded] = placeQueens(board,col+1); if succeeded return else % Backtrack. Undo the placement of the queen % on this column and keep looking. board(row,col) = 0; end end end

% If we get here, we have checked every row on this column % without succeeding. Return with failure. succeeded = false;

While I was generating the animations for boards of different sizes, two things caught my eye. First, there's no solution for a 6x6 board that has a queen in any corner. Look at the animated solution above and see if you can convince yourself about that.

Second, I noticed that a solution for the 5x5 board can be found without any backtracking. Check out the animation:

Next time, I'll write about the graphics programming and GIF and video file writing I did to do the algorithm visualization.

Get
the MATLAB code

Published with MATLAB® R2017a

Reader Uri Pe'er found an old post from 2011 about shortest paths in binary images, and he asked this follow-up question (paraphrased):... read more >>

]]>Reader Uri Pe'er found an old post from 2011 about shortest paths in binary images, and he asked this follow-up question (paraphrased):

*How can I determine if there is a path from one row to another in a binary image? And can I constrain the direction of movement along the path?*

Last week I showed how to tackle this problem using `imfill`. Today I'm going to use a `graph`.

In December 2015 I wrote about how to create graphs based on the connectivity of image pixels, and I submitted some functions to the File Exchange to make it easy.

One of those functions is `binaryImageGraph`. Let's try it using one of the same images from my last post.

```
bw1 = imread('http://blogs.mathworks.com/steve/files/path-shapes-1.png');
imshow(bw1)
```

g = binaryImageGraph(bw1)

g = graph with properties: Edges: [180444×2 table] Nodes: [45880×3 table]

Each foreground pixel in the binary image corresponds to a graph node, so the total number of graph nodes is the same as the number of foreground pixels.

nnz(bw1)

ans = 45880

The function `plotImageGraph` (from the same File Exchange submission) can be used to visualize the image graph.

plotImageGraph(g)

Zoom in to see the graph details in a small portion of the image.

axis([275 300 125 140])

Now let's figure out whether there are any graph paths connecting a row-100 pixels to a row-300 pixel. Note that `binaryImageGraph` stores the original foreground pixel locations in the graph's node table for purposes just like this. Here's what the node table looks like.

head(g.Nodes)

ans = 8×3 table x y PixelIndex ___ __ __________ 121 76 57676 121 77 57677 121 78 57678 121 79 57679 121 80 57680 121 81 57681 121 82 57682 121 83 57683

row_100_nodes = find(g.Nodes.y == 100); row_300_nodes = find(g.Nodes.y == 300);

Next, use `conncomp` to find the connected components of the graph.

bins = conncomp(g);

For one of the graph nodes, `bins(node)` tells you which connected component that node belongs to.

unique(bins(row_100_nodes))

ans = 1

So the foreground pixels on row 100 all belong to the same connected component.

unique(bins(row_300_nodes))

ans = 1 2

The foreground pixels on row 300 belong to two different connected components. Since one of those the components is the same as for the row 100 pixels, that means there is path from row 100 to row 300.

Just like `imfill`, `binaryImageGraph` can take a custom connectivity matrix. We can use that to constrain the allowed paths. Let's try that with one of the other images I used last time.

```
bw3 = imread('http://blogs.mathworks.com/steve/files/path-shapes-3.png');
imshow(bw3)
```

conn = [ 1 0 0 0 1 0 0 0 1 ]; g = binaryImageGraph(bw3,conn); plotImageGraph(g) axis([275 300 125 140])

Is there a constrained path from row 180 to row 250?

row_180_nodes = find(g.Nodes.y == 180); row_250_nodes = find(g.Nodes.y == 250); bins = conncomp(g); ~isempty(intersect(bins(row_180_nodes),bins(row_250_nodes)))

ans = logical 0

No.

But what if we reverse the path constraint?

conn = [ 0 0 1 0 1 0 1 0 0 ]; g = binaryImageGraph(bw3,conn); plotImageGraph(g) axis([275 300 125 140])

row_180_nodes = find(g.Nodes.y == 180); row_250_nodes = find(g.Nodes.y == 250); bins = conncomp(g); ~isempty(intersect(bins(row_180_nodes),bins(row_250_nodes)))

ans = logical 1

Yes, rows 180 and 250 are connected via a paths constrained to lie along the upper-right to lower-left diagonal.

Get
the MATLAB code

Published with MATLAB® R2017a

Reader Uri Pe'er found an old post from 2011 about shortest paths in binary images, and he asked this follow-up question (paraphrased):... read more >>

]]>Reader Uri Pe'er found an old post from 2011 about shortest paths in binary images, and he asked this follow-up question (paraphrased):

*How can I determine if there is a path from one row to another in a binary image? And can I constrain the direction of movement along the path?*

That sounds fun. Let's give it a try.

I can think of several different approaches. Today I'm going to tackle it using `imfill`.

Here's a sample image I just drew up.

```
bw1 = imread('http://blogs.mathworks.com/steve/files/path-shapes-1.png');
imshow(bw1)
```

Let's use `imfill` to see if there is a path from any of the foreground pixels on row 100 to any of the foreground pixels on row 300. First, find the foreground pixels on those rows.

fg_cols_row_100 = find(bw1(100,:)); fg_cols_row_300 = find(bw1(300,:)); hold on plot(fg_cols_row_100,100,'b*') plot(fg_cols_row_300,300,'b*') hold off

Now let's perform an `imfill` operation starting from all the foreground pixels on row 100. Because `imfill` fills in background pixels, we need to complement the image first.

```
bw1a = ~bw1;
imshow(bw1a)
title('Complemented image')
```

Now perform the fill in the complemented image.

```
seed_pixels = [100*ones(numel(fg_cols_row_100),1) fg_cols_row_100'];
bw1b = imfill(bw1a,seed_pixels);
imshow(bw1b)
title('Filled image')
```

Which pixels changed as a result of the fill?

```
bw1c = bw1b & ~bw1a;
imshow(bw1c)
title('Pixels changed by the fill operation')
```

If there are any foreground pixels on row 300 of `bw1c`, then that means row 100 is connected to row 300.

connected = any(bw1c(300,:))

connected = logical 1

Let's try the steps above on another image.

bw2 = imread('http://blogs.mathworks.com/steve/files/path-shapes-2.png'); imshow(bw2) title('Image 2')

fg_cols_row_100 = find(bw1(100,:)); fg_cols_row_300 = find(bw1(300,:)); bw2a = ~bw2; seed_pixels = [100*ones(numel(fg_cols_row_100),1) fg_cols_row_100']; bw2b = imfill(bw2a,seed_pixels); bw2c = bw2b & ~bw2a; connected = any(bw2c(300,:))

connected = logical 0

So we have computed that there is no path connecting rows 100 and 300 in the second image.

OK, what about the second part of Uri's question? How do we constrain the path? For example, how do we constrain the problem to allow movement from the upper-left-to-lower-right diagonal?

We can do that by specifying a *custom connectivity* when we call `imfill`. Here's another sample image to test the concept.

bw3 = imread('http://blogs.mathworks.com/steve/files/path-shapes-3.png'); imshow(bw3) fg_cols_row_100 = find(bw3(100,:)); fg_cols_row_300 = find(bw3(300,:)); hold on plot(fg_cols_row_100,100,'b*') plot(fg_cols_row_300,300,'b*') hold off

You can specify a custom connectivity using a 3x3 matrix. Here is one that only allows movement along the upper-left-to-lower-right diagonal.

conn = [ 1 0 0 0 1 0 0 0 1 ];

Now use that connectivity matrix with `imfill`.

bw3a = ~bw3; seed_pixels = [100*ones(numel(fg_cols_row_100),1) fg_cols_row_100']; bw3b = imfill(bw3a,seed_pixels,conn); bw3c = bw3b & ~bw3a; imshow(bw3c)

The odd-looking shape above is the set of pixels reachable from foreground pixels on row 100 traveling only along the specified diagonal. With this constrained path, rows 100 and 300 are not connected.

connected = any(bw3c(300,:))

connected = logical 0

But what if we use the opposite diagonal for our path constraint?

conn = [ 0 0 1 0 1 0 1 0 0 ]; bw3b = imfill(bw3a,seed_pixels,conn); bw3c = bw3b & ~bw3a; imshow(bw3c)

You can see that rows 100 and 300 are still not connected. Let's try rows 180 and 250.

imshow(bw3) fg_cols_row_180 = find(bw3(180,:)); fg_cols_row_250 = find(bw3(250,:)); hold on plot(fg_cols_row_180,180,'b*') plot(fg_cols_row_250,250,'b*') hold off

seed_pixels = [180*ones(numel(fg_cols_row_180),1) fg_cols_row_180']; bw3b = imfill(bw3a,seed_pixels,conn); bw3c = bw3b & ~bw3a; imshow(bw3c)

And rows 180 and 250 are connected.

connected = any(bw3c(250,:))

connected = logical 1

Next time I plan to take another look at this problem using MATLAB `graph` operations.

Get
the MATLAB code

Published with MATLAB® R2017a

One of my favorite aspects of MATLAB for image processing is how I can use a binary image to index directly into a grayscale image. The technique is called logical indexing, and I'm going to show you how it works today.... read more >>

]]>One of my favorite aspects of MATLAB for image processing is how I can use a binary image to index directly into a grayscale image. The technique is called *logical indexing*, and I'm going to show you how it works today.

*Note: this is an update of a post I originally wrote in 2008*.

Let me start with a small example. (As regular readers know, I like to use magic squares for small matrix examples.)

A = magic(5)

A = 17 24 1 8 15 23 5 7 14 16 4 6 13 20 22 10 12 19 21 3 11 18 25 2 9

Every MATLAB user is familiar with ordinary matrix indexing notation.

A(2,3)

ans = 7

`A(2,3)` extracts the 2nd row, 3rd column of the matrix `A`. You can extract more than one row and column at the same time:

A(2:4, 3:5)

ans = 7 14 16 13 20 22 19 21 3

When an indexing expression appears on the left-hand side of the equals sign, that's an assigment. You are changing one or more of the values of the variable on the left-hand side.

A(5,5) = 100

A = 17 24 1 8 15 23 5 7 14 16 4 6 13 20 22 10 12 19 21 3 11 18 25 2 100

Here is a frequently-asked MATLAB question: *How do I replace all the NaNs in my matrix B with 0s?*

An experienced MATLAB user will immediately answer:

B(isnan(B)) = 0;

For example:

B = rand(3,3); B(2, 2:3) = NaN

B = 0.2217 0.3188 0.0855 0.1174 NaN NaN 0.2967 0.5079 0.8010

Replace the NaNs with zeros:

B(isnan(B)) = 0

B = 0.2217 0.3188 0.0855 0.1174 0 0 0.2967 0.5079 0.8010

The expression `B(isnan(B))` is an example of *logical indexing*. Logical indexing is a compact and expressive notation that's very useful for many image processing operations.

Let's talk about the basic rules of logical indexing, and then we'll reexamine the expression `B(isnan(B))`.

If `C` and `D` are matrices, then `C(D)` is a logical indexing expression if `D` is a *logical* matrix.

*Logical* is one of the fundamental data types for MATLAB arrays. Relational operators, such as `==` or `>`, produce logical arrays automatically.

C = hilb(4)

C = 1.0000 0.5000 0.3333 0.2500 0.5000 0.3333 0.2500 0.2000 0.3333 0.2500 0.2000 0.1667 0.2500 0.2000 0.1667 0.1429

D = C > 0.4

D = 4×4 logical array 1 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0

If we use `D` as an index into `C` with the expression `C(D)`, then we will extract all the values of `C` corresponding to nonzero values of `D` and returns them as a column vector. It is equivalent to `C(find(D))`.

C(D)

ans = 1.0000 0.5000 0.5000

Now we know enough to break down the `B(isnan(B))` expression to see how it works.

B = rand(3,3); B(2, 2:3) = NaN; nan_locations = isnan(B)

nan_locations = 3×3 logical array 0 0 0 0 1 1 0 0 0

B(nan_locations)

ans = NaN NaN

B(nan_locations) = 0

B = 0.0292 0.4886 0.4588 0.9289 0 0 0.7303 0.2373 0.5468

Functions in the Image Processing Toolbox, as well as the MATLAB functions `imread` and `imwrite`, follow the convention that logical matrices are treated as binary (black and white) images. For example, when you read a 1-bit image file using `imread`, it returns a logical matrix:

bw = imread('text.png'); whos bw

Name Size Bytes Class Attributes bw 256x256 65536 logical

This convention, together with logical indexing, makes it very convenient and expressive to use binary images as pixel masks for extracting or operating on sets of pixels.

Here's an example showing how to use logical indexing to compute the histogram of a subset of image pixels. Specifically, given a grayscale image and a binary segmentation, compute the histogram of just the foreground pixels in the image.

Here's our original image:

```
I = imread('rice.png');
imshow(I)
```

Here's a segmentation result (computed and saved earlier), represented as a binary image:

```
url = 'http://blogs.mathworks.com/images/steve/192/rice_bw.png';
bw = imread(url);
imshow(bw)
```

Now use the segmentation result as a logical index into the original image to extract the foreground pixel values.

```
foreground_pixels = I(bw);
whos foreground_pixels
```

Name Size Bytes Class Attributes foreground_pixels 17597x1 17597 uint8

Finally, compute the histogram of the foreground pixels.

figure imhist(foreground_pixels)

As another example, you could complement the binary image to compute something based on the background pixels.

imhist(I(~bw))

*PS. I expect that this will be the last blog post that I write using R2016b. Keep an eye on the downloads page!*

Get
the MATLAB code

Published with MATLAB® R2016b

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