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

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