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.
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.
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.
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.
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:
>> 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:
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.
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:
the basic idea of finding shortest paths by adding two distance transforms together (part 1)
using bwdistgeodesic to find shortest paths subject to constraint (part 5)
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.
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:
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.
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.
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.
(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.)
"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.
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.
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.
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.
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:
Compute the distance transform for just the upper left block of foreground pixels.
Compute the distance transform for just the lower right block of foreground pixels.
Add the two distance transforms together.
Round to lower precision.
The pixels in the regional minimum of the result lie along one or more of the shortest paths from one block to the other.
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)
Earlier this month I started exploring the question of computing the shortest path from one object in a binary image to another. I described
the following basic algorithm:
Compute the distance transform for just the upper left block of foreground pixels.
Compute the distance transform for just the lower right block of foreground pixels.
Add the two distance transforms together.
The pixels in the regional minimum of the result lie along one or more of the shortest paths from one block to the other.
This algorithm only works for path-based approximations to the Euclidean distance transform; it doesn't work for the Euclidean
distance transform itself.
The Image Processing Toolbox function supports three path-based approximations to the Euclidean distance transform: 'cityblock',
'chessboard', and 'quasi-euclidean'.
Today I want to compare these three and discuss ambiguities associated with the idea of "shortest path" on an image.
Step 1 of the algorithm is to compute the distance transform for the first object. I'll use the 'cityblock' distance transform
approximation. In this approximation, only horizontal and vertical path segments are allowed. Diagonal segments are not allowed.
The distance between a pixel and any of its N, E, S, or W neighbors is 1.
The smallest value of D, 12, is the minimum path length (for paths consisting only of horizontal and vertical segments) from one object to the other.
Step 4 is to find the set of pixels in the regional minimum of D. This set represents the shortest path.
paths = imregionalmin(D);
To help visualize the path I'll use imoverlay from the MATLAB Central File Exchange. The code below will overlay the pixels on the shortest path in gray.
P = false(size(bw));
P = imoverlay(P, paths, [.5 .5 .5]);
P = imoverlay(P, bw, [1 1 1]);
imshow(P, 'InitialMagnification', 'fit')
Uh oh, why do we have a big rectangular block of gray instead of a single path? The answer is that there isn't a unique shortest
path. There are many paths you could travel to get from one object to the other that all have the same path length, 12. Here
are just a view of them:
subplot(2,2,1)
imshow(P, 'InitialMagnification', 'fit')
x = [2 6 6];
y = [2 2 10];
hold on
plot(x, y, 'g', 'LineWidth', 2)
hold off
subplot(2,2,2)
imshow(P, 'InitialMagnification', 'fit')
x = [2 2 6];
y = [2 10 10];
hold on
plot(x, y, 'g', 'LineWidth', 2);
hold off
subplot(2,2,3)
imshow(P, 'InitialMagnification', 'fit')
x = [2 3 3 4 4 5 5 6 6 6 6 6 6];
y = [2 2 3 3 4 4 5 5 6 7 8 9 10];
hold on
plot(x, y, 'g', 'LineWidth', 2)
hold off
subplot(2,2,4)
imshow(P, 'InitialMagnification', 'fit')
x = [2 2 2 2 2 2 3 3 4 4 5 5 6];
y = [2 3 4 5 6 7 7 8 8 9 9 10 10];
hold on
plot(x, y, 'g', 'LineWidth', 2)
hold off
Next, let's look at the 'chessboard' distance. In this distance transform approximation, a path can consist of horizontal,
vertical, and diagonal segments, and the path length between a pixel and any of its neighbors (N, NE, E, SE, S, SW, W, and
NW) is 1.0.
Now the set of gray pixels is not rectangular, but it still encompasses multiple possible shortest paths. And, surprisingly,
it seems to include some pixels on which the path actually moves temporarily away from its destination.
Here are several possible paths, all of length 8 (according to the 'chessboard' distance transform).
subplot(2,2,1)
imshow(P, 'InitialMagnification', 'fit')
x = [2 3 4 5 6 6 6 6 6];
y = [2 3 4 5 6 7 8 9 10];
hold on
plot(x, y, 'g', 'LineWidth', 2)
hold off
subplot(2,2,2)
imshow(P, 'InitialMagnification', 'fit')
x = [2 2 2 2 2 3 4 5 6];
y = [2 3 4 5 6 7 8 9 10];
hold on
plot(x, y, 'g', 'LineWidth', 2);
hold off
subplot(2,2,3)
imshow(P, 'InitialMagnification', 'fit')
x = [2 1 2 1 2 3 4 5 6];
y = [2 3 4 5 6 7 8 9 10];
hold on
plot(x, y, 'g', 'LineWidth', 2)
hold off
subplot(2,2,4)
imshow(P, 'InitialMagnification', 'fit')
x = [2 3 4 3 2 3 4 5 6];
y = [2 3 4 5 6 7 8 9 10];
hold on
plot(x, y, 'g', 'LineWidth', 2)
hold off
Finally, let's look at the 'quasi-euclidean' distance transform. In this variation, paths from one pixel to another can consist
of horizontal, vertical, and diagonal segments. The distance from a pixel to one of its horizontal or vertical neighbors is
1, while the distance to one of its diagonal neighbors is sqrt(2).
Here's the computation again, this time specifying 'quasi-euclidean':
I am trying to find the length and location of the shortest path between objects. It seems like some application of bwdist
should be able to do this. Any advice? Example:
1 0 1
0 1 0
0 0 0
0 1 0
1 0 1
Should yield
0 0 0
0 0 0
0 p 0
0 0 0
0 0 0
With p = path length = 1
I agree, some application of bwdist should do it. Actually, there are several other useful functions that can be applied to this task, including imregionalmin and bwdistgeodesic, a new function that just shipped in R2011b.
I've been tinkering with this problem off and on for a few days, and I think it might take a few posts to explore some of
the interesting ideas and issues posed by this problem. Today I'll start with the basics of using bwdist and imregionalmin to identify a set of shortest paths from one object to another.
Let's use this simple test image to start our exploration.
We want to find the shortest path (or paths!) from the upper left block of foreground pixels to the lower right block. The
key algorithmic idea is this:
Compute the distance transform for just the upper left block of foreground pixels.
Compute the distance transform for just the lower right block of foreground pixels.
Add the two distance transforms together.
The pixels in the regional minimum of the result lie along one or more of the shortest paths from one block to the other.
(For reference, see P. Soille, Morphological Image Analysis: Principles and Applications, 2nd edition, Springer, 2003, section
7.3.4, "Application to minimal path detection," pp. 233-234.)
But we can't use the normal Euclidean distance transform because it doesn't correspond to a path traversal from one pixel
center to another. Rather, we have to use one of the distance transform approximations supported by bwdist: 'chessboard', 'cityblock', and 'quasi-euclidean'.
You can see that the smallest value of D_sum is about 7.8284, and there's a set of pixels in-between the two objects that all have that value. The function imregionalmin identifies them all:
In the quasi-euclidean approximation to the distance transform, the distance from a pixel to one of its horizontal or vertical
neighbors is 1, and the distance from a pixel to one of its diagonal neighbors is sqrt(2). The result above says that the
shortest such path between the two objects is 7.8284. Furthermore, it says that there is more than one shortest path.
I plan to explore this topic more fully this month. When using the quasi-euclidean distance transform, for example, there
are floating-point round-off issues that have to be dealt with when you try to apply the technique above to larger images.
Another issue is the effect of using different path-based distance transform approximations, such as 'cityblock' or 'chessboard'. Yet another issue is how to identify specific paths from the family of paths identified by imregionalmin. Along the way I'll experiment with different ways to visualize the results.
Finally, I'll explore how to use a new function, bwdistgeodesic, to find shortest paths when there are constraints on where the paths can go.
All the posts in this series
the basic idea of finding shortest paths by adding two distance transforms together (part 1)
Steve Eddins is a software development manager in the MATLAB and image processing areas at MathWorks. Steve coauthored Digital Image Processing Using MATLAB. He writes here about image processing concepts, algorithm implementations, and MATLAB.
Recent Comments