I started this blog at the beginning of 2006. At first, I left comments open and only filtered out obvious spam. Soon, though, I found that I couldn't handle the volume of comments. Also, most comments were off-topic, which I thought detracted from the main content of the posts. So I began moderating all comments, and I only allowed through comments that were relevant to the post's topic.

And that's the way it's been for the past five years or so. But lately the time required for handling and responding to comments has become a problem for me. It is taking away time and energy that I would prefer to put into posting new content. Because of this I have occasionally considered turning off comments completely, but I've been reluctant to do that because I regularly learn valuable things from reader comments, and also because readers sometimes let me know about corrections I should make.

So starting this week I'm trying a compromise. I'm closing comments on any post older than 60 days. I believe that most of the interesting discussions and needed corrections come in well before that time.

Keep in mind that this blog isn't a general help forum or an alternative channel for product support. I encourage everyone seeking specific help with a MATLAB concept or with getting your code running to use resources such as MATLAB Answers or the MATLAB newsgroup.

For product support questions, another starting point is the Support section of mathworks.com. There's really a lot of good information there.

I also get a pretty high volume of direct e-mail because my e-mail address is relatively easy to find. I do try to answer specific questions about my MATLAB Central File Exchange submissions that come in through my Author page there. Otherwise, I mostly can't answer the direct e-mail that I receive.

Thanks for reading ... and for commenting!

I was hired to take over development of the Image Processing Toolbox, which had just shipped version 1.0. From the summary, I see that I spent a lot of time that first month studying the toolbox, looking at the way product demos and GUIs were programmed (back then, anyway), and writing a new product demo (on quadtree decomposition).

A fair amount of time that first month was spent getting up and running on my computer, an HP workstation running some flavor of Unix. I read a few months worth of internal newsgroup postings in order to learn more about my new colleagues and company. I drafted a new section on JPEG image compression for a marketing document called MATLAB Expo, and I searched for a cool-looking sample image to use in an ad. I attended a meeting of the MathWorks Social Mission Team.

I studied specs for MATLAB 5.0, which we thought would ship soon but actually took three more years.

And I learned the company's bug tracking system in order to report my very first bug found.

Hey ... speaking of new jobs, we just completed our annual planning process, and in 2012 we plan to hire bunches of people at MathWorks!

So please, take a look at the Careers section of our web site to see if an opening interests you. And submit your resume now! Managers are going crazy trying to fill the openings.

This week I came across some files I wrote about 16 years ago to compute Hilbert curves. A Hilbert curve is a type of fractal curve; here is a sample:

I can't remember why I was working on this. Possibly I was anticipating that 16 years in the future, during an unusually mild New England winter, I would be looking for a blog topic.

Anyway, there are several interesting ways to code up a Hilbert curve generator. My old code for generating the Hilbert curve followed the J. G. Griffiths, "Table-driven Algorithms for Generating Space-Filling Curves," Computer-Aided Design, v.17, pp. 37-41, 1985. Here's the basic procedure:

First, initialize four curves, \( A_0 \), \( B_0 \), \( C_0 \), and \( D_0 \), to be empty (no points).

Then, build up the Hilbert curve iteratively as follows:

\[ A_{n+1} = [B_n, N, A_n, E, A_n, S, C_n] \]

\[ B_{n+1} = [A_n, E, B_n, N, B_n, W, D_n] \]

\[ C_{n+1} = [D_n, W, C_n, S, C_n, E, A_n] \]

\[ D_{n+1} = [C_n, S, D_n, W, D_n, N, B_n] \]

where \( N \) represents a unit step up, \( E \) is a unit step to the right, \( S \), is a unit step down, and \( W \) is a unit step left.

One way to code this procedure is to incrementally build up a set of vectors that define the step from one point on the path to the next, and then to use a cumulative summation at the end to turn the steps into x-y coordinates. Here's how you might do it.

A = zeros(0,2);
B = zeros(0,2);
C = zeros(0,2);
D = zeros(0,2);

north = [ 0  1];
east  = [ 1  0];
south = [ 0 -1];
west  = [-1  0];

order = 3;
for n = 1:order
  AA = [B ; north ; A ; east  ; A ; south ; C];
  BB = [A ; east  ; B ; north ; B ; west  ; D];
  CC = [D ; west  ; C ; south ; C ; east  ; A];
  DD = [C ; south ; D ; west  ; D ; north ; B];

  A = AA;
  B = BB;
  C = CC;
  D = DD;
end

A = [0 0; cumsum(A)];

plot(A(:,1), A(:,2), 'clipping', 'off')
axis equal, axis off

I was curious to see what might be on the MATLAB Central File Exchange, so I searched for "hilbert curve" and found several interesting contributions. There are a couple of 3-D Hilbert curve generators, and several different ways of coding up a 2-D Hilbert curve generator. I was particularly interested in the Fractal Curves contribution by Jonas Lundgren.

Jonas shows a much more compact implementation of the ideas above using complex arithmetic. It looks like this:

a = 1 + 1i;
b = 1 - 1i;

% Generate point sequence
z = 0;
order = 5;
for k = 1:order
    w = 1i*conj(z);
    z = [w-a; z-b; z+a; b-w]/2;
end

plot(z, 'clipping', 'off')
axis equal, axis off

Jonas' contribution includes several other curve generators. Here's one called "dragon".

z = dragon(12);
plot(z), axis equal

Here's a zoom-in view of a portion of the dragon.

axis([-.1 0.2 0.4 0.7])

What do you think? Does anyone know of an image processing application for shape-filling curves? Please leave a comment.

Get the MATLAB code

Published with MATLAB® 7.13

Much of the information I posted in this blog years ago is still useful today. Image processing theory hasn't been completely been overturned since then, and I'm still talking about MATLAB after all. For the benefit of readers who have joined the party more recently, I will occasionally recap the useful bits that I posted about five years ago.

In August, September, and October of 2006 I finally finished (mostly) a long series of posts covering spatial transforms. I posed a reader challenge to come up with a particularly interesting and creative custom spatial transformation. My favorite submission was this one:

I told the story of how we guessed wrong in the early 1990s about which of the several variants of the discrete cosine transform we should use in the Image Processing Toolbox.

I poked fun at advertising techniques for digital cameras.

I taught about several image processing concepts, including Hough transforms, separable convolution, and nonflat grayscale dilation and erosion.

And I finally ended ten years of public silence about the MATLAB default image.

Be sure to check out the blog archives for more oldies but goodies.

Get the MATLAB code

Published with MATLAB® 7.13

Blog reader Mike posed the following question recently:

If you have a bunch of point locations (for example, object centroids), how you make a binary image containing just those points?

For example, consider this image:

bw = imread('text.png');
imshow(bw, 'InitialMagnification', 200)

How can we make an image like this, where the dots are located at the centroids of the objects?

Solving this problem is a nice application of linear indexing, something I wrote about in this blog a long time ago. Let's see how it can work for us here.

First, let's find the centroids using regionprops:

s = regionprops(bw, 'Centroid');

s is a struct array. Since we just asked for one measurement, the centroid, each element of s is a struct containing just one field, 'Centroid'.

s(1)

ans = 
    Centroid: [11 13.5000]

s(2)

ans = 
    Centroid: [7.6829 38.1707]

The length of s is the number of objects in the image.

num_objects = length(s)

num_objects =
    88

Next, we gather all the individual centroid locations into x and y vectors. To accomplish this I use the comma-separated list syntax for struct arrays.

centroids = cat(1, s.Centroid);
x = centroids(:,1);
y = centroids(:,2);

If the comma-separated list syntax makes your brain hurt, you can use a loop instead:

  centroids = zeros(length(s), 2);
  for k = 1:length(s)
      centroids(k,:) = s(k).Centroid;
  end

Now let's round the centroid locations to get row and column subscripts.

r = round(y);
c = round(x);

Here's where linear indexing comes into play. In order to assign to a bunch of scattered locations like this, you want to use a single subscript. That's what we call linear indexing. You can use the function sub2ind to convert a set of subscripts to linear indices.

ind = sub2ind(size(bw), r, c);

And finally we can use the linear indices to assign a value to a bunch of image pixel locations all at once.

bw2 = false(size(bw));
bw2(ind) = true;

imshow(bw2, 'InitialMagnification', 200)

See my 08-Feb-2008 blog post for more about linear indexing.

Get the MATLAB code

Published with MATLAB® 7.13

Blog reader Asadullah posted the following question last week on my old post about batch processing:

I am trying to process some images by following the MATLAB demo. After getting the names of files when I try to see any of the files then it gives the error. The detail is as follows:

 >> fileFolder = fullfile(matlabroot,'work','original');
 >> dirOutput = dir(fullfile(fileFolder,'AT3_1m4_*.jpg'));
 >> fileNames = {dirOutput.name}'

 fileNames =

     'AT3_1m4_001.jpg'
     'AT3_1m4_002.jpg'
     'AT3_1m4_003.jpg'
     'AT3_1m4_004.jpg'
     'AT3_1m4_005.jpg'
     'AT3_1m4_006.jpg'
     'AT3_1m4_007.jpg'

     [...]

     'AT3_1m4_116.jpg'
     'AT3_1m4_117.jpg'
     'AT3_1m4_118.jpg'
     'AT3_1m4_119.jpg'
     'AT3_1m4_120.jpg'

 >> I=imread(fileNames{1});
 ??? Error using ==> imread
 File "AT3_1m4_001.jpg" does not exist.

I could not understand where is the problem? Please help me to solve this.

Several people have asked me this same question, so I thought I should post the answer so everyone can see it instead of replying in a comment.

First, though, I want to strongly recommend against placing your own work files (data, code, etc.) underneath the MATLAB program folder, as the Asadullah's code shows in this line:

   fileFolder = fullfile(matlabroot,'work','original');

The function matlabroot returns the location of the MATLAB program folder. On a typical Windows computer, for example, this location might be something like:

   C:\Program Files\MATLAB\R2011b

Windows doesn't consider this to be an area where a user will store their own files, and most backup programs will not normally back up files that are stored here. Instead, store your own work elsewhere. On Windows, a good choice is somewhere under "My Documents".

OK, now back to the question at hand. The problem is that the string being passed to imread contains only the filename, such as 'AT3_1m4_116.jpg', and not the full folder location. That's not enough information for imread to be able to find the file.

You have to add the full folder location to the image filename before passing it to imread. In Asadullah's example, the corrected code would look something like this:

   I = imread(fullfile(fileFolder,fileNames{1}));

Hope that helps! May your days be filled with happily processing large numbers of images.

Get the MATLAB code

Published with MATLAB® 7.13

For the past several weeks I've been writing about shortest-path problems in image processing: finding the shortest path between two points in an image, with and without constraints.

Applications included the practical (path finding in a skeleton image) and fun (maze solving).

Along the way, I've described:

I hope you found the series interesting and useful. If you have particular applications for these techniques in your own work, I invite you to share them with us by posting a comment.

Get the MATLAB code

Published with MATLAB® 7.13

In this post in the Exploring shortest paths series, I make things a little more complicated (and interesting) by adding constraints to the shortest path problem. Last time, I showed this example of finding the shortest paths between the "T" and the sideways "s" in this image:

url = 'http://blogs.mathworks.com/images/steve/2011/two-letters-cropped.png';
bw = imread(url);

And this was the result:

What if we complicate the problem by adding other objects to the image and looking for the shortest path that avoids these other objects?

For example, what if we want to find a shortest path between the "T" and the sideways "s" without going through any of other letters between them in the image below?

url = 'http://blogs.mathworks.com/images/steve/2011/text-cropped.png';
text = imread(url);
imshow(text, 'InitialMagnification', 'fit')

We can do this by using bwdistgeodesic instead of bwdist in our algorithm. bwdistgeodesic is a new function in the R2011b release of the Image Processing Toolbox. In addition to taking a binary image input that bwdist takes, it takes another argument that specifies which pixels the paths are allowed to traverse.

Let's use our text image to make a mask of allowed path pixels.

mask = ~text | bw;
imshow(mask, 'InitialMagnification', 'fit')

Next, we make our two binary object images as before.

L = bwlabel(bw);
bw1 = (L == 1);
bw2 = (L == 2);

Next, instead of using bwdist, we use bwdistgeodesic.

D1 = bwdistgeodesic(mask, bw1, 'quasi-euclidean');
D2 = bwdistgeodesic(mask, bw2, 'quasi-euclidean');

D = D1 + D2;
D = round(D * 8) / 8;

D(isnan(D)) = inf;
paths = imregionalmin(D);

paths_thinned_many = bwmorph(paths, 'thin', inf);
P = false(size(bw));
P = imoverlay(P, ~mask, [1 0 0]);
P = imoverlay(P, paths, [.5 .5 .5]);
P = imoverlay(P, paths_thinned_many, [1 1 0]);
P = imoverlay(P, bw, [1 1 1]);
imshow(P, 'InitialMagnification', 'fit')

In the image above, the white characters are the starting and ending points for the path. The path is not allowed to pass through the red characters. The gray pixels are all the pixels on any shortest path, and the yellow pixels lie along one particular shortest path.

Let's have a little fun by using this technique to solve a maze. Here's a test maze (created by The Maze Generator).

url = 'http://blogs.mathworks.com/images/steve/2011/maze-51x51.png';
I = imread(url);
imshow(I)

Make an image of walls by thresholding the maze image. Dilate the walls (make them a little thicker) in order to keep the solution path a little bit away from them.

walls = imdilate(I < 128, ones(7,7));
imshow(walls)

The function bwgeodesic will take seed locations (as column and row coordinates) in addition to binary images. Below I specify the seed locations as the two entry points into the maze.

D1 = bwdistgeodesic(~walls, 1, 497, 'quasi-euclidean');
D2 = bwdistgeodesic(~walls, 533, 517, 'quasi-euclidean');

The rest of the procedure is the same as above.

D = D1 + D2;
D = round(D * 8) / 8;

D(isnan(D)) = inf;
paths = imregionalmin(D);

solution_path = bwmorph(paths, 'thin', inf);
thick_solution_path = imdilate(solution_path, ones(3,3));
P = imoverlay(I, thick_solution_path, [1 0 0]);
imshow(P, 'InitialMagnification', 'fit')

(If you're interested in other maze-solving techniques using image processing, take a look at the File Exchange contribution maze_solution by Image Analyst.)

In a comment on part 1 of this series, Sepideh asked this question:

"I have a binary image which contains one pixel width skeleton of an image. I have got the skeleton using the bwmorph function. As you may know basically in this binary image, I have ones on the skeleton and zeros anywhere else. how can I find the shortest ‘quasi-euclidean’ distance between two nonzero elements which are on the skeleton? The path should be on the skeleton itself too."

The techniques described above can be applied to this path-along-a-skeleton problem. Let's try an example.

bw = imread('circles.png');
imshow(bw, 'InitialMagnification', 200)

skel = bwmorph(bw, 'skel', Inf);
imshow(skel, 'InitialMagnification', 200)

Let's find the skeleton-constrained path distance between these two points.

r1 = 38;
c1 = 53;

r2 = 208;
c2 = 147;

hold on
plot(c1, r1, 'g*', 'MarkerSize', 15)
plot(c2, r2, 'g*', 'MarkerSize', 15)
hold off

D1 = bwdistgeodesic(skel, c1, r1, 'quasi-euclidean');
D2 = bwdistgeodesic(skel, c2, r2, 'quasi-euclidean');

D = D1 + D2;
D = round(D * 8) / 8;

D(isnan(D)) = inf;
skeleton_path = imregionalmin(D);
P = imoverlay(skel, imdilate(skeleton_path, ones(3,3)), [1 0 0]);
imshow(P, 'InitialMagnification', 200)
hold on
plot(c1, r1, 'g*', 'MarkerSize', 15)
plot(c2, r2, 'g*', 'MarkerSize', 15)
hold off

Finally, the skeleton path length is given by the value of D at any point along the path.

path_length = D(skeleton_path);
path_length = path_length(1)

path_length =

  239.2500

All the posts in this series

• the basic idea of finding shortest paths by adding two distance transforms together (part 1)
• the nonuniqueness of the shortest paths (part 2)
• handling floating-point round-off effects (part 3)
• using thinning to pick out a single path (part 4)
• using bwdistgeodesic to find shortest paths subject to constraint (part 5)

Get the MATLAB code

Published with MATLAB® 7.13

In my previous posts on Exploring shortest paths (November 1, November 26, and December 3), I have noted several times that there isn't a unique shortest path (in general) between one object and another in a binary image. To illustrate, here's a recap of the 'quasi-euclidean' example from last time.

bw = logical([ ...
    0 0 0 0 0 0 0
    0 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 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 0 0 0 1 0
    0 0 0 0 0 0 0  ]);
L = bwlabel(bw);
bw1 = (L == 1);
bw2 = (L == 2);

D1 = bwdist(bw1, 'quasi-euclidean');
D2 = bwdist(bw2, 'quasi-euclidean');
D = D1 + D2;
D = round(D * 32) / 32;

paths = imregionalmin(D);
P = false(size(bw));
P = imoverlay(P, paths, [.5 .5 .5]);
P = imoverlay(P, bw, [1 1 1]);
imshow(P, 'InitialMagnification', 'fit')

The gray pixels shown above all belong to one or more of the shortest paths between the two objects. So how do we pick one path?

The first idea I had was to thin the paths "blob" using the 'thin' option of bwmorph.

paths_thinned_once = bwmorph(paths, 'thin');
P = false(size(bw));
P = imoverlay(P, paths, [.5 .5 .5]);
P = imoverlay(P, paths_thinned_once, [1 1 0]);
P = imoverlay(P, bw, [1 1 1]);
imshow(P, 'InitialMagnification', 'fit')

That seems promising. What if we keep thinning until convergence? You do that by passing inf as the repetition value to bwmorph.

paths_thinned_many = bwmorph(paths, 'thin', inf);
P = false(size(bw));
P = imoverlay(P, paths, [.5 .5 .5]);
P = imoverlay(P, paths_thinned_many, [1 1 0]);
P = imoverlay(P, bw, [1 1 1]);
imshow(P, 'InitialMagnification', 'fit')

That looks great.

Let's try the whole sequence with a longer path.

url = 'http://blogs.mathworks.com/images/steve/2011/two-letters-cropped.png';
bw = imread(url);
imshow(bw, 'InitialMagnification', 'fit')

L = bwlabel(bw);
bw1 = (L == 1);
bw2 = (L == 2);

D1 = bwdist(bw1, 'quasi-euclidean');
D2 = bwdist(bw2, 'quasi-euclidean');
D = D1 + D2;
D = round(D * 32) / 32;

paths = imregionalmin(D);

paths_thinned_many = bwmorph(paths, 'thin', inf);
P = false(size(bw));
P = imoverlay(P, paths, [.5 .5 .5]);
P = imoverlay(P, paths_thinned_many, [1 1 0]);
P = imoverlay(P, bw, [1 1 1]);
imshow(P, 'InitialMagnification', 'fit')

As before, the pixels belonging to any shortest path are shown in gray. The pixels in yellow belong to one particular path.

Next time I'll explain how to find shortest paths subject to constraints.

All the posts in this series

• the basic idea of finding shortest paths by adding two distance transforms together (part 1)
• the nonuniqueness of the shortest paths (part 2)
• handling floating-point round-off effects (part 3)
• using thinning to pick out a single path (part 4)
• using bwdistgeodesic to find shortest paths subject to constraint (part 5)

Get the MATLAB code

Published with MATLAB® 7.13

In part 2 of Exploring shortest paths, I noted a problem with using the 'quasi-euclidean' distance transform to find the shortest paths between two objects in a binary image. Specifically, our algorithm resulted in a big, unfilled gap between the two objects.

bw = logical([ ...
    0 0 0 0 0 0 0
    0 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 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 0 0 0 1 0
    0 0 0 0 0 0 0  ]);
L = bwlabel(bw);
bw1 = (L == 1);
bw2 = (L == 2);

D1 = bwdist(bw1, 'quasi-euclidean');
D2 = bwdist(bw2, 'quasi-euclidean');
D = D1 + D2;

paths = imregionalmin(D);
P = false(size(bw));
P = imoverlay(P, paths, [.5 .5 .5]);
P = imoverlay(P, bw, [1 1 1]);
imshow(P, 'InitialMagnification', 'fit')

The pixels marked in gray should be the set of pixels that lie along a shortest path from the two objects in the original image. You can plainly, see, however, that there's a large gap.

To see why, let's examine more closely the values of D.

D

D =

   12.4853   11.6569   11.6569   12.2426   12.8284   13.4142   14.8284
   11.0711    9.6569   10.2426   10.8284   11.4142   12.0000   13.4142
   10.4853    9.6569    9.6569   10.2426   10.8284   11.4142   12.8284
   10.4853    9.6569    9.6569    9.6569   10.2426   10.8284   12.2426
   10.4853    9.6569    9.6569    9.6569    9.6569   10.2426   11.6569
   11.0711    9.6569    9.6569    9.6569    9.6569    9.6569   11.0711
   11.6569   10.2426    9.6569    9.6569    9.6569    9.6569   10.4853
   12.2426   10.8284   10.2426    9.6569    9.6569    9.6569   10.4853
   12.8284   11.4142   10.8284   10.2426    9.6569    9.6569   10.4853
   13.4142   12.0000   11.4142   10.8284   10.2426    9.6569   11.0711
   14.8284   13.4142   12.8284   12.2426   11.6569   11.6569   12.4853

There appears to be an unbroken set of pixels between the two objects with value 9.6569. Actually, it turns out that the pixels do not all have the same value. For example, D(3,3) appear to have the same value D(4,4) but are actually slightly different.

D(3,3)

ans =

    9.6569

D(4,4)

ans =

    9.6569

D(3,3) - D(4,4)

ans =

 -9.5367e-007

The difference between the two is the corresponding single-precision relative floating-point precision:

eps(D(3,3))

ans =

  9.5367e-007

This floating-point round-off error comes into play because of the way that multiples of sqrt(2) are being added in different orders.

So we need to adjust our algorithm to account for floating-point round-off error. One way to do that is to round the distance transform values to some lower precision. For example, the code below rounds the distance transform to be a multiple of (1/32).

D = round(D*32)/32;

Now D(3,3) and D(4,4) are equal:

D(3,3) - D(4,4)

ans =

     0

And imregionalmin works as expected to extract the set of pixels belonging to shortest paths. (I'm using imoverlay from the MATLAB Central File Exchange.)

paths = imregionalmin(D);
P = false(size(bw));
P = imoverlay(P, paths, [.5 .5 .5]);
P = imoverlay(P, bw, [1 1 1]);
imshow(P, 'InitialMagnification', 'fit')

Here's our revised algorithm, including the new rounding step:

Next time I'll look into how to pick a particular path among the many shortest-path choices.

All the posts in this series

• the basic idea of finding shortest paths by adding two distance transforms together (part 1)
• the nonuniqueness of the shortest paths (part 2)
• handling floating-point round-off effects (part 3)
• using thinning to pick out a single path (part 4)
• using bwdistgeodesic to find shortest paths subject to constraint (part 5)

Get the MATLAB code

Published with MATLAB® 7.13