This is the first of a few blog posts about object measurements based on a concept called the Feret diameter, sometimes called the caliper diameter. The diagram below illustrates the concept. Place the object to be measured inside the jaws of a caliper, with the caliper oriented at a specified... read more >>

]]>This is the first of a few blog posts about object measurements based on a concept called the *Feret diameter*, sometimes called the *caliper diameter*. The diagram below illustrates the concept. Place the object to be measured inside the jaws of a caliper, with the caliper oriented at a specified angle. Close the jaws tightly on the object while maintaining that angle. The distance between the jaws is the Feret diameter at that angle.

Eventually, we'll examine several related measurements (and algorithms for computing them), like these.

Today, though, I'm just going to play around with a very small binary image. (Note: the `pixelgrid` function is a utility that I wrote. It's at the bottom of this post.)

```
bw = [ ...
0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 1 0
0 0 0 0 0 0 1 1 1 1 0
0 0 0 0 0 1 1 1 1 0 0
0 0 0 1 1 1 1 1 1 0 0
0 0 1 1 1 0 0 0 0 0 0
0 0 1 1 1 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 ];
imshow(bw)
pixelgrid
```

I want to brute-force my to finding the two corners of the white blob that are farthest apart. Let's start by finding all the pixel corners. (I'll be using implicit expansion here, which works in R2016b or later.)

[y,x] = find(bw); hold on plot(x,y,'xb','MarkerSize',5) corners = [x y] + cat(3,[.5 .5],[.5 -.5],[-.5 .5],[-.5 -.5]); corners = permute(corners,[1 3 2]); corners = reshape(corners,[],2); corners = unique(corners,'rows'); plot(corners(:,1),corners(:,2),'sr','MarkerSize',5) hold off

Well, let's get away from pure brute force just a little. We really only need to look at the points that are part of the *convex hull* of all the corner points.

hold on k = convhull(corners); hull_corners = corners(k,:); plot(hull_corners(:,1),hull_corners(:,2),'r','LineWidth',3) plot(hull_corners(:,1),hull_corners(:,2),'ro','MarkerSize',10,'MarkerFaceColor','r') hold off

Now let's find which two points on the convex hull are farthest apart.

dx = hull_corners(:,1) - hull_corners(:,1)'; dy = hull_corners(:,2) - hull_corners(:,2)'; pairwise_dist = hypot(dx,dy); [max_dist,j] = max(pairwise_dist(:)); [k1,k2] = ind2sub(size(pairwise_dist),j); point1 = hull_corners(k1,:); point2 = hull_corners(k2,:); hold on plot([point1(1) point2(1)],[point1(2) point2(2)],'-db','LineWidth',3,'MarkerFaceColor','b') hold off

And the maximum distance is ...

max_dist

max_dist = 10

(I thought I had made a mistake when I first saw this number come out to be an exact integer. But no, the answer is really 10.)

That's a pretty good start. I don't know exactly where I'll go with this topic next time. I'm making it up as I go.

function hout = pixelgrid(h) %pixelgrid Superimpose a grid of pixel edges on an image % pixelgrid superimposes a grid of pixel edges on the image in the % current axes. % % pixelgrid(h_ax) superimposes the grid on the image in the specified % axes. % % pixelgrid(h_im) superimposes the grid on the specified image. % % group = pixelgrid(___) returns an hggroup object that contains the two % lines that are used to draw the grid. One line is thick and darker, and % the other is thin and lighter. Using two contrasting lines in this way % guarantees that the grid will be visible regardless of pixel colors. % % EXAMPLES % % Superimpose pixel grid on color image. Zoom in. % % rgb = imread('peppers.png'); % imshow(rgb) % pixelgrid % axis([440 455 240 250]) % % Change the colors and line widths of the pixel grid. % % rgb = imread('peppers.png'); % imshow(rgb) % h = pixelgrid; % axis([440 455 240 250]) % h.Children(1).Color = [178 223 138]/255; % h.Children(1).LineWidth = 2; % h.Children(2).Color = [31 120 180]/255; % h.Children(2).LineWidth = 4; % % LIMITATIONS % % This function is intended for use when looking at a zoomed-in image % region with a relatively small number of rows and columns. If you use % this function on a typical, full-size image without zooming in, the % image will not be visible under the grid lines. % % This function attempts to determine if MATLAB is running with a high % DPI display. If so, then it displays the pixel grid using thinner % lines. However, the method for detecting a high DPI display on works on % a Mac when Java AWT is available. % Steve Eddins % Copyright The MathWorks, Inc. 2017 if nargin < 1 him = findobj(gca,'type','image'); elseif strcmp(h.Type,'axes') him = findobj(h,'type','image'); elseif strcmp(h.Type,'image') him = h; else error('Invalid graphics object.') end if isempty(him) error('Image not found.'); end hax = ancestor(him,'axes'); xdata = him.XData; ydata = him.YData; [M,N,~] = size(him.CData); if M > 1 pixel_height = diff(ydata) / (M-1); else % Special case. Assume unit height. pixel_height = 1; end if N > 1 pixel_width = diff(xdata) / (N-1); else % Special case. Assume unit width. pixel_width = 1; end y_top = ydata(1) - (pixel_height/2); y_bottom = ydata(2) + (pixel_height/2); y = linspace(y_top, y_bottom, M+1); x_left = xdata(1) - (pixel_width/2); x_right = xdata(2) + (pixel_width/2); x = linspace(x_left, x_right, N+1); % Construct xv1 and yv1 to draw all the vertical line segments. Separate % the line segments by NaN to avoid drawing diagonal line segments from the % bottom of one line to the top of the next line over. xv1 = NaN(1,3*numel(x)); xv1(1:3:end) = x; xv1(2:3:end) = x; yv1 = repmat([y(1) y(end) NaN], 1, numel(x)); % Construct xv2 and yv2 to draw all the horizontal line segments. yv2 = NaN(1,3*numel(y)); yv2(1:3:end) = y; yv2(2:3:end) = y; xv2 = repmat([x(1) x(end) NaN], 1, numel(y)); % Put all the vertices together so that they can be drawn with a single % call to line(). xv = [xv1(:) ; xv2(:)]; yv = [yv1(:) ; yv2(:)]; hh = hggroup(hax); dark_gray = [0.3 0.3 0.3]; light_gray = [0.8 0.8 0.8]; if isHighDPI bottom_line_width = 1.5; top_line_width = 0.5; else bottom_line_width = 3; top_line_width = 1; end % When creating the lines, use AlignVertexCenters to avoid antialias % effects that would cause some lines in the grid to appear brighter than % others. line(... 'Parent',hh,... 'XData',xv,... 'YData',yv,... 'LineWidth',bottom_line_width,... 'Color',dark_gray,... 'LineStyle','-',... 'AlignVertexCenters','on'); line(... 'Parent',hh,... 'XData',xv,... 'YData',yv,... 'LineWidth',top_line_width,... 'Color',light_gray,... 'LineStyle','-',... 'AlignVertexCenters','on'); % Only return an output if requested. if nargout > 0 hout = hh; end end function tf = isHighDPI % Returns true if running with a high DPI default display. % % LIMITATION: Only works on a Mac with Java AWT available. In other % situations, this function always returns false. tf = false; if ismac if usejava('awt') sf = java.awt.GraphicsEnvironment.getLocalGraphicsEnvironment().getDefaultScreenDevice.getScaleFactor(); tf = sf > 1; end end end

Get
the MATLAB code

Published with MATLAB® R2017b

I happen to know the author, and he's an OK fellow.

]]>Before I dive into today's topic, I want to invite you to come back to the MATLAB Central Blogs page in a couple of weeks. Look for something new to appear related to deep learning.... read more >>

]]>*Before I dive into today's topic, I want to invite you to come back to the MATLAB Central Blogs page in a couple of weeks. Look for something new to appear related to deep learning*.

Earlier today, I was looking at recent image processing questions on MATLAB Answers, and I came across this question:

"Using regionprops to change intensity of blobs"

I thought to myself, "I bet I've written about something like before." It turns out that I have -- but it was more than ten years ago! It is ever-so-slightly possible that a couple of you have not been reading this blog that long, so I thought I would update the post for you today.

The function `regionprops` is very useful for measuring the properties of shapes in a binary image. There are documentation examples and product demos showing how to do this, and I've shown this function in action several times in this blog.

But sometimes I get questions about how to process pixel values in the "original" gray scale image. In other words, suppose your process is something like this:

- Segment gray scale image to get a binary image of objects
- Do something with the original gray scale pixel values corresponding to each blob in the binary image

This post is about how to do that last step. Let's start by creating some simple sample images to work with.

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

bw = imbinarize(I); bw = imfill(bw, 'holes'); imshow(bw) title('Thresholded, with holes filled')

The key to processing the pixels in `I` using the blobs in `bw` is to use `regionprops` to get the `'PixelIdxList'` for each blob.

```
s = regionprops(L, 'PixelIdxList')
```

s = 10×1 struct array with fields: PixelIdxList

s(1)

ans = struct with fields: PixelIdxList: [2651×1 double]

`s` is a struct array. `s(k).PixelIdxList` is a vector of the linear indices of the k-th region. You can use these values to index directly into `I`. For example, here is the mean pixel value of the 3rd coin:

idx = s(3).PixelIdxList; mean(I(idx))

ans = 172.8769

You can also replace object pixel values in `I` by indexing on the left side of an assignment statement, like this:

```
% Turn the 4th coin white and the 5th coin black:
J = I;
J(s(4).PixelIdxList) = 255;
J(s(5).PixelIdxList) = 0;
imshow(J)
```

Here's a loop that replaces the pixel values for each coin by the mean pixel value for that coin:

F = I; for k = 1:numel(s) idx = s(k).PixelIdxList; F(idx) = mean(I(idx)); end imshow(F)

I'll finish with a question that Mat asked in a blog comment: "How can I find the location of the maximum pixel value in each labeled region?" Here's how to find the locations and then plot them on top of the original image. (In addition to `'PixelIdxList'`, I'll also ask `regionprops` for `'PixelList'`, which gives the x and y coordinates of each pixel in the region.)

s = regionprops(L, 'PixelIdxList', 'PixelList'); imshow(I) hold on for k = 1:numel(s) idx = s(k).PixelIdxList; pixels_in_region_k = I(idx); [max_value, max_idx] = max(pixels_in_region_k); max_location = s(k).PixelList(max_idx, :); plot(max_location(1),max_location(2),'o','MarkerSize',10,... 'MarkerFaceColor','w','MarkerEdgeColor','k') end hold off

Get
the MATLAB code

Published with MATLAB® R2017a

Today I want to tell you how and why I made these images:... read more >>

]]>Today I want to tell you how and why I made these images:

After the MATLAB R2014b release, I wrote several blog posts (Part 1, Part 2, Part 3, and Part 4) about the new default colormap, parula, that was introduced in that release.

Sometime later, I came across some material by Peter Kovesi about designing perceptually uniform colormaps (or colourmaps, as Peter writes it).

I was particularly intrigued by a test image that Peter designed for the purpose of visually evaluating the perceptual characteristics of a colormap. Here is the test image:

The image is constructed by superimposing a sinusoid on a linear ramp, with the amplitude of the sinusoid getting smaller as you move away from the top row. Here are three cross-sections of the image: row 1, row 64, and row 128.

url = 'http://peterkovesi.com/projects/colourmaps/cm_MATLAB_gray256.png'; I = im2double(imread(url)); subplot(3,1,1) plot(I(1,:)) axis([1 512 0 1]) title('Row 1') subplot(3,1,2) plot(I(64,:)) axis([1 512 0 1]) title('Row 64') subplot(3,1,3) plot(I(128,:)) axis([1 512 0 1]) title('Row 128')

Here is the basic code for making this test image. I'm going to vary Kovesi's image slightly. I'll add an extra half-cycle of the sinusoid so that it reaches a peak at the right side of the image, and I'll add a full-range linear ramp section at the bottom. (If you look closely at the cross-section curves above, you'll see that the linear ramp goes from 5% to 95% of the range.)

First, compute the ramp.

num_cycles = 64.5; pixels_per_cycle = 8; A = 0.05; width = pixels_per_cycle * num_cycles + 1; height = round((width - 1) / 4); ramp = linspace(A, 1-A, width);

Next, compute the sinusoid.

k = 0:(width-1); x = -A*cos((2*pi/pixels_per_cycle) * k);

Now, vary the amplitude of the sinusoid with the square of the distance from the bottom of the image.

q = 0:(height-1); y = ((height - q) / (height - 1)).^2; I1 = (y') .* x;

Superimpose the sinusoid on the ramp.

I = I1 + ramp;

Finally, add a full-range linear ramp section to the bottom of the image.

```
I = [I ; repmat(linspace(0,1,width), round(height/4), 1)];
clf
imshow(I)
title('Colormap test image')
```

Last week, I posted Colormap Test Image to the File Exchange. It contains the function `colormapTestImage`, which does all this for you.

I = colormapTestImage;

The function has another syntax, too. If you pass it the name of a colormap, it will display the test image using that colormap. For example, here is the test image with the old MATLAB default colormap, `jet`.

```
colormapTestImage('jet')
```

This test image illustrates why we replaced `jet` as the default MATLAB colormap. I have annotated the image below to show some of the issues.

Now compare with the new default colormap, `parula`.

```
colormapTestImage('parula')
```

I think that illustrates what we were trying to achieve with `parula`: perceptual fidelity to the data.

Since I'm talking about `parula`, I'll finish by mentioning that we need some very subtle tweaks to `parula` in the R2017b release. So you can compare, I'll show you the original version that shipped with R2014b.

```
colormapTestImage('parula_original')
```

Readers, can you tell what is different? Let us know in the comments.

Get
the MATLAB code

Published with MATLAB® R2017a

Wow.... read more >>

]]>I saw a presentation last month that mentioned a user request to have the ability to customize regionprops. That is, a user wanted to be able to add their own measurement to regionprops.... read more >>

]]>I saw a presentation last month that mentioned a user request to have the ability to customize `regionprops`. That is, a user wanted to be able to add their own measurement to `regionprops`.

Today, I'll show you how to do this yourself.

First, here's a brief recap on what `regionprops` does. The function computes measurements of image regions. Some of these measurements are based purely on a region's shape, while others incorporate pixel values within the regions. Here's an example using the coins.png sample image.

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

Let's convert this image to binary, using adaptive thresholding, filling holes, and removing small "noise" pixels.

bw = imbinarize(I,'adaptive'); bw = imfill(bw,'holes'); bw = bwareafilt(bw,[100 Inf]); imshow(bw)

You can count the "blobs" (object) yourself; there are 10 of them.

The simplest `regionprops` call, `regionprops(bw)` computes the `Area`, `Centroid`, and `BoundingBox` for each object.

s = regionprops(bw)

s = 10×1 struct array with fields: Area Centroid BoundingBox

But I don't think this is the best way to call `regionprops` anymore. You can now tell `regionprops` to return the results as a table.

```
t = regionprops('table',bw)
```

t = 10×3 table Area Centroid BoundingBox ____ ________________ ____________ 2635 37.133 106.85 [1x4 double] 1846 56.131 49.693 [1x4 double] 2672 96.199 146.05 [1x4 double] 1839 109.97 84.848 [1x4 double] 2744 120.37 208.73 [1x4 double] 2520 148.57 34.404 [1x4 double] 2589 174.83 120.01 [1x4 double] 2518 216.81 70.649 [1x4 double] 1857 236.03 173.36 [1x4 double] 1829 265.96 102.64 [1x4 double]

The table form is a lot more convenient for many tasks. For today's topic, one especially nice thing thing about tables is how easy it is to add your own variables to the table.

To illustrate, let's add a measurement that I've seen called `Roundness`. One definition for roundness is:

$R = \frac{4A}{\pi L^2}$

where $A$ is the object area and $L$ is the major axis length of the best-fit ellipse for the object. Here's how to compute roundness and add it directly to the measurements returned by `regionprops`.

First, note that both `Area` and `MajorAxisLength` are supported by `regionprops`, so let's start with those.

t = regionprops('table',bw,'Area','MajorAxisLength')

t = 10×2 table Area MajorAxisLength ____ _______________ 2635 60.08 1846 50.178 2672 59.792 1839 49.674 2744 60.374 2520 58.08 2589 58.676 2518 58.162 1857 49.77 1829 49.564

You access table variables using dot notation, like `t.Area`. Similarly, you can create a new table variable using dot notation and assignment, like `t.MyVariable = ...`. So adding `Roundness` to the table returned by `regionprops` is this simple.

t.Roundness = 4 * t.Area ./ (pi * t.MajorAxisLength.^2)

t = 10×3 table Area MajorAxisLength Roundness ____ _______________ _________ 2635 60.08 0.92945 1846 50.178 0.93352 2672 59.792 0.9516 1839 49.674 0.94893 2744 60.374 0.9585 2520 58.08 0.95118 2589 58.676 0.95745 2518 58.162 0.94772 1857 49.77 0.95453 1829 49.564 0.94798

Let's try this computation with an image containing objects that are not quite as round.

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

bw2 = imbinarize(I2,'adaptive'); bw2 = imfill(bw2,'holes'); bw2 = bwareafilt(bw2,[100 Inf]); imshow(bw2)

t2 = regionprops('table',bw2,'Area','MajorAxisLength'); t2.Roundness = 4 * t2.Area ./ (pi * t2.MajorAxisLength.^2); head(t2)

ans = 8×3 table Area MajorAxisLength Roundness ____ _______________ _________ 138 23.594 0.31562 120 18.152 0.4637 169 28.123 0.27207 157 23.793 0.3531 284 43.757 0.18885 200 26.259 0.36929 141 21.647 0.38311 177 29.087 0.26636

I'm a big fan of the (relatively) new `histogram` function in MATLAB, so let's use it to compare our roundness numbers. I will follow the advice given in the `histogram` reference page for normalizing multiple histograms so that they can be more easily compared. I'll set the y-axis limits to [0 1], which is appropriate for probability normalization, and I'll set the x-axis limits to [0 1], which is the range for `Roundness`.

h1 = histogram(t.Roundness); hold on h2 = histogram(t2.Roundness); hold off h1.Normalization = 'probability'; h2.Normalization = 'probability'; h1.BinWidth = 0.02; h2.BinWidth = 0.02; xlim([0 1]); ylim([0 1]); title('Histogram of roundness (probability normalization)') legend('coins','rice')

There you have it. You can add your own object measurements to the output of `regionprops`. It's especially easy if you tell `regionprops` to return a table.

I'll leave you with this question, dear reader: Are there measurements you would like us to add to `regionprops`? I am aware of an enhancement request for the Feret diameter. What else would you like to see?

Get
the MATLAB code

Published with MATLAB® R2017a

A MATLAB user recently contacted tech support with a question about the output of `fft2`. The user had a function, q(x,y), evaluated on an (x,y) grid. The grid evaluation produced a matrix that the user passed to `fft2`. The user's question: what are the spatial frequencies associated with the elements of the output from `fft2`?

The question was forwarded to the MATLAB Math team. Chris Turnes, a fellow Georgia Tech grad, answered it. Here is Chris' answer, lightly edited.

*To map discrete Fourier coefficients into samples of the corresponding Discrete Time Fourier Transform, we need to know about the sampling pattern of the spatial signal. Since the customer presumably has the vectors x and y that she used to query the continuous function q(x,y), this ought to be easy to determine.*

*To determine the sampling frequencies, we use the spacing of the sampling grid. The sampling frequencies in each dimension are fsx = 1/(x(2)-x(1)) and fsy = 1/(y(2)-y(1)) if the grid is uniform in each dimension.*

*Next, the normalized DFT frequencies of a (shifted) N-point DFT (where N here is even for convenience) are fa = ((-N/2):(N/2-1))/N. Therefore, the corresponding DTFT frequency grid would be composed of the frequencies fcx = fa*fsx and fcy = fa*fsy. That would be the correct way to do the mapping between FFT coefficients and the corresponding continuous spatial frequency.*

Thanks, Chris.

Get
the MATLAB code

Published with MATLAB® R2017a

A few weeks ago, I wrote about a solver algorithm for the Eight Queens Problem. The post included diagrams like this one.... read more >>

]]>A few weeks ago, I wrote about a solver algorithm for the Eight Queens Problem. The post included diagrams like this one.

Today I want to show you how I made that diagram in MATLAB. First, let's talk about the board. Back in 2011, I wrote about a variety of ways to make a checkerboard (or chessboard) pattern. In that post, I played games with `repmat`, powers of -1, `floor`, and `round`. It got a little crazy.

I'm still fond of using integer powers of -1, like this.

$(-1)^{(n_1 + n_2)}$

In 2011, I used `ndgrid` to make a matrix of integer powers, $n_1 + n_2$.

[n1,n2] = ndgrid(0:4); (-1).^(n1 + n2)

ans = 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1

Since the introduction of implicit expansion in R2016b, though, I no longer need to use `ndgrid` to explicitly form the matrix of powers.

n = 0:4; (-1).^(n + n')

ans = 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1

In my eight queens code, here's how I got the shades of gray for a chessboard. (There are lots of ways to do this.)

N = 8; dark_square_color = .7; light_square_color = .9; color_range = light_square_color - dark_square_color; f = (-1).^((1:N)' + (1:N)); f = (f + 1) * color_range / 2; f = f + dark_square_color; f(1:4,1:4)

ans = 0.9000 0.7000 0.9000 0.7000 0.7000 0.9000 0.7000 0.9000 0.9000 0.7000 0.9000 0.7000 0.7000 0.9000 0.7000 0.9000

Next, I needed to replicate each element of `f` to make an image with larger squares. Do you know how to replicate elements of a matrix? Some experienced MATLAB users would do it using the `kron` function. Now, however, you can just use the new `repelem` function. The reference page tells you when this function was introduced.

I'm going to make my squares 60 pixels wide, and then I'll turn it into a truecolor image so that the pixel colors are independent of the figure's colormap.

```
pixel_size = 60;
f = repelem(f,pixel_size,pixel_size);
f = repmat(f,1,1,3);
h = imshow(f,'InitialMagnification',100);
```

Now I want to show you a coordinate system trick. If you turn on the axes display, you can see the image pixel coordinates. (In Image Processing Toolbox terminology, these are called *intrinsic coordinates*.)

```
axis on
```

When I get to the part about displaying the queens, however, I would like to be able to place them based on the coordinates of each chessboard square, independent of the number of pixels per square. I can do that by manipulating the `XData` and `YData` properties of the image. These properties assign spatial coordinates to the left/right and top/bottom edges of the image.

h.XData = [0.5 N+0.5]; h.YData = [0.5 N+0.5]; axis([0.5 N+0.5 0.5 N+0.5])

Now let's put some queens on the board. This turns out to be pretty easy because you can just do it with a `text` object. Since R2014b, you've been able to draw text in MATLAB graphics using Unicode characters. And the Unicode character set includes chess symbols! The black queen symbol is at code point 9819 (decimal).

queen = char(9819)

queen = '♛'

Let's put a queen on the square in the second row, third column. (That's $x=3$, $y=2$.)

hq = text(3,2,queen);

Oops, our queen is pretty small, and she also appears to be off-center. Let's fix that.

```
hq.FontSize = 40;
hq.HorizontalAlignment = 'center';
```

That's better.

I have one more little coding trick to show you. When I wrote my eight queens animation code, I wanted to have an easy to show and remove a queen at any square on the board. So, I made an matrix of text objects, one at each square. Then I could just index into the matrix of text objects and turn the `Visible` property on and off. Here's how. (First, let me delete the text object I just made.)

delete(hq) for r = 1:N for c = 1:N queens(r,c) = text(c,r,char(9819),... 'HorizontalAlignment','center',... 'Visible','off',... 'FontSize',40); end end

OK, let's turn on the queen displays on a couple of squares.

queens(5,3).Visible = 'on'; queens(2,4).Visible = 'on';

For a final bit of fun, the function `set` lets you modify properties of a whole array of graphics objects at once.

Turn them all on!

set(queens,'Visible','on')

(The board above, by the way, is NOT a solution to the Eight Queens Problem. In case you were wondering.)

Get
the MATLAB code

Published with MATLAB® R2017a

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