Steve on Image Processing

Today's post is about this picture, which is from the paper "Classroom Examples of Robustness Problems in Geometric Computations," by L. Kettner, K. Mehlhorn, S. Pion, S. Schirra, and C. Yap, Computational Geometry, vol. 40, n. 1, May 2008, pp. 61-78.

Yes, I know, I was supposed to do another follow-up on my Cody problem of eliminating unnecessary polygon vertices. But I really am following up on the Cody problem, if a bit indirectly.

I was looking at solution 108897 by Cris, and I paused over these lines:

v1 = p2-p1;
v2 = p3-p2;
if ( v1(1)*v2(2) == v1(2)*v2(1) ) && any(v1.*v2>0)

I recognized the expression ( v1(1)*v2(2) == v1(2)*v2(1) ) as a way to see if the point p2 is on the line segment formed by p1 and p3. It reminded me of Loren's 06-Jun-2008 blog post on methods for determining collinearity, and it also reminded me of the paper above.

I decided to see if I could reproduce the figure in MATLAB.

The figure concerns the behavior of floating-point methods for computing the planar orientation predicate. Given three points $p$, $q$, and $r$ in a plane, does the point $r$ lie to the left of the line through $p$ and $q$, to the right of that line, or on the line? One way to compute this is via the equation

$$O(p,q,r) = \mbox{sign}((q_x - p_x)(r_y - p_y) - (q_y - p_y)(r_x - p_x))$$

The figure in the paper represents an experiment using $q=(12,12)$, $r=(24,24)$, and $256^2$ different locations for the point $p$, all of which are extremely close to the point $p=(0.5,0.5)$. More specifically, the figure is a visualization of $O((p_x+xu, p_y+yu),q,r)$ for $0 \leq x,y \leq 255$ and $u=2^{-53}$. This value of $u$ is the increment between adjacent double-precision floating-point numbers in the interval $0.5 < x \leq 1$.

Here's a straightforward implementation of the orientation equation above.

function out = floatOrient(p,q,r)

out = sign((q(1) - p(1))*(r(2) - p(2)) - ...
    (q(2) - p(2))*(r(1) - p(1)));

Let's see what we get when we compute the orientation for all those different point locations $p'$ that are extremely close to $(0.5,0.5)$.

[m,n] = meshgrid(0:255);
Px = 0.5 + n*eps(0.5);
Py = 0.5 + m*eps(0.5);

q = [12 12];
r = [24 24];

out = zeros(size(Px));
for u = 1:256
    for v = 1:256
        out(u,v) = floatOrient([Px(u,v), Py(u,v)], q, r);
    end
end

(Yes, I know ... floatOrient and the calling code above can all be vectorized.)

I'll show in yellow the locations $p'$ where $r$ has been determined to be on the line containing $p'$ and $q$. I'll use blue and red to show where $r$ has been determined to lie on one side or the other of the line.

image(out + 2);
axis equal;
set(gca,'ydir','normal')
axis off
colormap([0 0 1; 1 1 0; 1 0 0])
hold on
plot([.5 256.5],[.5 256.6],'k','LineWidth',3)
hold off

How do we interpret this? First, I want to emphasize that this matrix represents the results for $p'$ that are all extremely close to $(0.5,0.5)$. The lower left corner corresponds exactly to the $(0.5,0.5)$, while the upper right corner corresponds to approximately $(0.500000000000028,0.500000000000028)$. Ideally, only the results lying along the black line should be yellow; these are the locations where $p'$, $q$, and $r$ are collinear. Our floating-point computation method is showing many other locations in yellow. Notice also that a few pairs of locations $p'$ are actually being computed reversed. A left turn has been determined to be a right turn, and vice versa.

In the comment thread on Loren's collinearity post, Tim Davis suggested using a rank test to determine collinearity. It looks something like this:

function out = rankCollinear(p,q,r)
out = rank([p-q ; q-r]) < 2;

Let's see how this test does on our collection of points.

for u = 1:256
    for v = 1:256
        out(u,v) = rankCollinear([Px(u,v), Py(u,v)], q, r);
    end
end

image(out + 1);
axis equal;
set(gca,'ydir','normal')
axis off
colormap([0 0 0; 1 1 0])
hold on
plot([.5 256.5],[.5 256.6],'k','LineWidth',3)
hold off

We only see yellow pixels, not red and blue, because rankCollinear only checks for collinearity, not the left or right turning of noncollinear points. The width of the section of yellow pixels depends on the numerical tolerance used by the rank function.

I am fascinated by this particular experiment and visualization. It reminds me once again that seemingly straightforward geometric ideas can be devilishly difficult to implement with a computer program.

Get the MATLAB code

Published with MATLAB® 7.14

Published with MATLAB® 7.14

A couple of weeks ago I posted about how one might want to simplify the output of bwboundaries in order to remove unneeded vertices in the output boundaries. That is, turn this polygon:

into this one:

At the same time, I created a new problem on Cody called "Eliminate unnecessary polygon vertices" and invited everyone to try it.

Well, in the last sixteen days or so, almost 300 solutions have been submitted for this problem!

I've been putting off writing a follow-up blog post because there's so much interesting to look at in the various solution attempts. But I decided today that I need to at least go ahead and recognize the efforts of everyone who participated in this effort.

(As an aside, the MATLAB Community team released updates to Cody just this week that make it a lot easier to see who has solved a problem and to search for all solutions by a particular player. That helped me a lot; thanks, team!)

By clicking on the "Solvers" button I can see the eight players who (as of this writing) have submitted correct solutions to this problem:

I especially appreciate the input from Sven, Richard, and Cris, who all helped "debug" my initial test suite in the first few hours.

Sven and Richard engaged in a very impressive tag team effort over the course of several days to find the smallest solution.

And I want to recognize Sagiv for making this his first solved Cody problem! (But keep going, Sagiv, try some more!)

When I was looking at various solutions, I found it helpful to create a function to help visualize them. You can download my visCody820Solution.m file from here. (You'll also need arrow.m by Erik Johnson.)

For example, here is the visualization of one of the test cases.

Next time I'll write about some of the interesting concepts and MATLAB techniques raised by this problem and its many solutions, including:

Get the MATLAB code

Published with MATLAB® 7.14

Published with MATLAB® 7.14

Today I'm combining my blog post with a new problem submission on Cody. The problem is related to the function bwboundaries in the Image Processing Toolbox. This function traces the boundaries of objects (and holes within objects), returning each boundary as a set of x-y vertices.

Here's an example from the documentation:

I = imread('rice.png');
BW = im2bw(I, graythresh(I));
[B,L] = bwboundaries(BW,'noholes');
imshow(label2rgb(L, @jet, [.5 .5 .5]))
hold on
for k = 1:length(B)
    boundary = B{k};
    plot(boundary(:,2), boundary(:,1), 'w', 'LineWidth', 2)
end
hold off

I have heard a few times from customers that they would like to eliminate "unnecessary vertices" in the output of bwboundaries. To illustrate, let's look at the boundary for a simple binary image containing a rectangle.

close(gcf)
bw = false(15,40);
bw(5:10,10:30) = true;
B = bwboundaries(bw);
boundary = B{1};
imshow(bw,'InitialMagnification', 'fit')
hold on
plot(boundary(:,2), boundary(:,1))
plot(boundary(:,2), boundary(:,1), '*')
hold off

Note that the boundary goes through the pixel centers, which is why a little bit of white appears outside the boundary line.

I can imagine that for some applications it would be nice to have the boundary polygon without as many vertices, like this:

x = [10 30 30 10 10];
y = [5 5 10 10 5];
imshow(bw,'InitialMagnification', 'fit')
hold on
plot(x, y)
plot(x, y, '*')
hold off

This idea has been on my potential blog topics list for a long time, and today I finally got around to thinking about it. I can think of several different approaches, and it's not immediately obvious to me which approach might be best (for whatever definition of best you believe in).

So I thought I would make a Cody problem out of it and see what creative ideas pop up.

I encourage you to give the problem a try. I'll summarize the results later.

Good luck!

UPDATE: In a week or so, I'll a pick the best solution (according to the well-known purple-seven multiphasic optimality criterion) and send out a prize (a hat or something similar) to the solver.

Get the MATLAB code

Published with MATLAB® 7.14

Published with MATLAB® 7.14

Last month I was experimenting with an algorithm in which I needed to construct a set of equidistant steps along a path consisting of line segments. (I'll call this a polyline path, or polyline for short.) Today I want to show you how to construct such a path in MATLAB. Along the way I'll show you the function improfile and give you a pointer to a database of sample range images that I found.

Contents

Polylines

Here's a simple example of a polyline. (Note that, because we're going to be talking about distances, I'm going to set the aspect ratio of my plots to be 1:1 by calling axis equal.)

x = [24 214 327 617];
y = [375 223 218 341];
plot(x,y)
axis equal

This has three connected line segments. In this example, each segment has a different length. I wanted to be able to find a certain number of points that were equally spaced along this path.

IMPROFILE

After pondering this problem for a minute or so, I remembered that the function improfile does exactly this computation internally. The function improfile is used to compute pixel-value profiles along a user-specified polyline.

Let me illustrate improfile using a range image. In a range image, the pixel values represent estimated distance from the imaging device. I found the USF (University of South Florida) Range Image Database online. Below is a range image I downloaded from the database. It was taken by the Vision Lab at USF using a "K2T structured light camera." (No, I don't know what that is.) The web site helpfully tells me that the image is stored in "rasterfile format," which is pretty dated but fortunately supported by imread.

url = 'http://blogs.mathworks.com/images/steve/2012/chessmen.range';
I = imread(url, 'ras');
imshow(I, 'InitialMagnification', 50)

Let's use improfile to compute a cross-section of range values along a particular path. First, though, let me superimpose the path on the image so you can where it is. (Although improfile has an interactive mode that lets you select the path using a mouse, I can't show that here.)

hold on
plot(x,y,'r','LineWidth',5)
hold off

Here's the cross-section of range values.

improfile(I,x,y)

You can see that improfile shows the cross-section as a 3-D plot. To see the cross-section as a normal 2-D line plot, call improfile with an output argument to get the cross-section values and then pass them to plot.

c = improfile(I,x,y);
plot(c)

Computing equidistant steps along the polyline

As part of its computation, improfile needs to compute equidistant steps along the polyline path. So how does that work?

The function first computes the cumulative distance from the beginning of the polyline to each vertex along the way, and then it uses a clever call to interp1 to compute the steps.

xy = [x' y'];
d = diff(xy,1)

d =

   190  -152
   113    -5
   290   123

dist_from_vertex_to_vertex = hypot(d(:,1), d(:,2))

dist_from_vertex_to_vertex =

  243.3187
  113.1106
  315.0063

cumulative_dist_along_path = [0;
    cumsum(dist_from_vertex_to_vertex,1)]

cumulative_dist_along_path =

         0
  243.3187
  356.4293
  671.4356

Now we can get our steps by constructing a call to interp1.

num_points = 20;
dist_steps = linspace(0, cumulative_dist_along_path(end), num_points);
points = interp1(cumulative_dist_along_path, xy, dist_steps)

points =

   24.0000  375.0000
   51.5949  352.9241
   79.1898  330.8482
  106.7847  308.7722
  134.3796  286.6963
  161.9745  264.6204
  189.5694  242.5445
  218.0483  222.8209
  253.3525  221.2587
  288.6567  219.6966
  323.9609  218.1345
  356.7328  230.6108
  389.2662  244.4095
  421.7996  258.2081
  454.3330  272.0068
  486.8664  285.8054
  519.3998  299.6041
  551.9332  313.4027
  584.4666  327.2014
  617.0000  341.0000

Plot the points, superimposing them on the original path, to see how we did.

plot(x,y)
hold on
plot(points(:,1), points(:,2), 'o')
hold off
axis equal

There you go. Now you know about improfile, you know how to compute a path along a polyline, and you know about an online database of sample range imagery.

I think we should call it a day.

Get the MATLAB code

Published with MATLAB® 7.14

Published with MATLAB® 7.14

On June 2 I showed you how to make a table with image colors and pixels appear when you publish your MATLAB scripts to HTML using the publish function. The result looks like this:

rgb = imread('peppers.png');
disp(im2html(rgb(88:92,200:204,:)))

Today I want to follow up with a few details about how this works.

First, I am taking advantage of the ability when publishing a MATLAB script to embed dynamically-generated HTML. To see how that works, let's start with some fixed HTML.

str = ['<html>' ...
    '<span style="color:white; background-color:black;">' ...
    'White text on black background' ...
    '</span>' ...
    '</html>'];
disp(str)

To make a simple dynamic example, let's make a random background color and then automatically assign the foreground color to be white or black, depending on how dark the background color is.

bg_color = rand(1,3)

bg_color =

    0.3404    0.5853    0.2238

bg_as_gray = rgb2gray(bg_color)

bg_as_gray =

    0.4709    0.4709    0.4709

if bg_as_gray(1) > 0.6
    fg_color = 'black';
else
    fg_color = 'white';
end

fg_color

fg_color =

white

Now make the HTML string again, this time inserting the dynamically generated colors. (First, though, convert the background color RGB values to be uint8 in the range [0,255].)

bg_color = im2uint8(bg_color);
str = sprintf(['<html>' ...
    '<span style="color:%s; background-color:rgb(%d,%d,%d);">' ...
    'Contrasting text on randomly colored background' ...
    '</span>' ...
    '</html>'], fg_color, bg_color(1), bg_color(2), bg_color(3));
disp(str)

So the im2html function simply produces an HTML table with dynamically generated table cell colors and pixel-value text strings. Here's an illustrative sample of the generated HTML:

<table class="image-table">
  <tr class="image-table-row">
    <!-- ROW: 1  COL: 1 -->
    <td class="dark-pixel-cell" style="background-color:#50324D; width:5.0em; height:5.0em">
      R:80<br \>G:50<br \>B:77
    </td>

Now we know the basics of dynamically generating some HTML and rendering it in your published output. How do we use that to generate our image pixel-color tables? It's more complicated than you might think.

First, remember that the MATLAB image display model is very flexible. The relationship between a pixel value and the actual color displayed on the screen can depend on these various factors:

(I wrote a whole series of blog posts about this topic way back when.)

Second, it can be complicated to figure out how to format the pixel-value text for all the different cases. Sometimes there is one value, sometimes three, sometimes both an index and a colormap value, sometimes the values integers, sometimes floating, ...

You get the idea.

Well, there's a tool in the Image Processing Toolbox called imagemodel that provides these services. We don't talk about it very much because we think it's not particularly useful to most toolbox users. However, if you are doing GUI programming involving images, it might come in handy.

imagemodel takes a Handle Graphics image handle and returns an object that you can query.

h = imshow('peppers.png');

model = imagemodel(h)

 
model =
 
IMAGEMODEL object accessing an image with these properties:

       ClassType: 'uint8'
    DisplayRange: []
     ImageHeight: 384
       ImageType: 'truecolor'
      ImageWidth: 512
    MinIntensity: []
    MaxIntensity: []

 

Here's how to get the RGB screen display color for a particular pixel:

getScreenPixelRGBValue(model, 5, 10)

ans =

    0.2588    0.1216    0.2510

I can get the RGB screen display color exactly the same way, even for a completely different kind of image. Below is a gray-scale image displayed with a black-to-white range of [-1,1].

theta = linspace(-pi,pi,200);
I = repmat(sin(theta),200,1);
h2 = imshow(I,[-1 1])

h2 =

  347.0363

model2 = imagemodel(h2);
getScreenPixelRGBValue(model2, 23, 47)

ans =

    0.0039    0.0039    0.0039

For getting pixel-value text strings, you have to ask the image model for a function handle that does the work. Then you can call the function handle. Here's how that would work for the sample image just above.

f = getPixelRegionFormatFcn(model2);

f(75, 128)

ans = 

    '0.76'

This is a fairly obscure corner of the Image Processing Toolbox, but I thought some of you might be interested to get a peek into how tools such as the Pixel Region Tool work.

Get the MATLAB code

Published with MATLAB® 7.14

Published with MATLAB® 7.14

All of us MathWorkers who hang out on MATLAB Central are very excited that MATLAB creator and MathWorks cofounder Cleve Moler has just become the latest MATLAB Central blogger.

I think I first started to get to know a little bit about Cleve from the Usenet newsgroup comp.soft-sys.matlab, which I started eagerly reading when it was formed at the beginning of 1993. Cleve posted "Hello from MathWorks" the first week, and he's been a very active participant ever since.

I met Cleve in person in September 1993 when I interviewed at MathWorks. I vaguely recall that he asked me if I knew anything about the use of wavelets and fractals for image compression, but mostly he wanted to show off the cool new secret stuff MathWorks was working on for the Symbolic Math Toolbox.

I mostly worked on the Image Processing Toolbox back then, but I sought Cleve's advice for one of my earliest contributions to MATLAB. I thought it might be beneficial for the filter2 function to automatically determine if the input filter is separable. (See my 04-Oct-2006 post on separable convolution.) Cleve explained to me how the rank function works and advised me on numerical issues and speed trade-offs. We decided to make the change, and it's been there ever since.

Cleve also helped me early on with an image processing algorithm. Someone suggested to me that we could fill in (or replace) pixels in a particular region by modeling it like a soap-film problem. That is, if you form wire in the shape of the pixel values (as heights) around the region to be filled, what shape would a soap film take inside the wire? I did a little searching (in books; the World Wide Web was still just a couple of years old) and found that I needed a solution where each pixel equals the average value of its north, east, south, and west neighbors, subject to boundary conditions based on the pixel values around the edge of the region. I programmed an iterative algorithm that worked, but it was quite slow.

When I asked Cleve for his advice, he said, "Isn't that a sparse system? Why don't you just use backslash?"

[head slap] Of course!

He was very nice about it, so I was only a little embarrassed.

I could go on and on about Cleve. However, he doesn't like it when people make a fuss about him, so he's probably already mad at me.

One more thought before I wrap this up. Anybody interested in how MATLAB came to be should watch Cleve's Origins of MATLAB video. He spins a great yarn.

Hey, Cleve - I've got two numbers whose average is three. What are they?

Get the MATLAB code

Published with MATLAB® 7.14

Published with MATLAB® 7.14

Today's post shows you how to make a table with image colors and pixels appear when you publish your MATLAB scripts to HTML using the publish function. The result looks like this:

rgb = imread('peppers.png');
disp(im2html(rgb(88:92,200:204,:)))

I was inspired to do something like this when I saw Printing Variables to HTML Tables in Published Code (by Ned) on the File Exchange a while back. It also produces an HTML table with colored cells and superimposed values. Here's a sample:

I was also thinking about the Pixel Region Tool in the Image Processing Toolbox. Here's a screen shot:

I wanted to go a bit further than Ned's original. I wanted to handle all the different kinds of image types (grayscale, truecolor, indexed with direct mapping, indexed with scaled mapping). I also wanted to replicate the feature of the Pixel Region Tool that automatically changed the color of the superimposed text depending on whether the underlying pixel was dark or light. (You can see that effect in the screen shot above.)

I've packaged all this in a function called im2html. You can download this function from the File Exchange.

Here are some examples showing how to use im2html with different types of images.

Display a table of values for a gray-scale image:

I = imread('pout.tif');
disp(im2html(I(125:134, 104:114)))

I = magic(10);
disp(im2html(I,[]))

Display a table of values from an indexed image:

[X,map] = imread('trees.tif');
disp(im2html(X(156:160,244:248),map))

You can also capture the output of im2html as a string, or write it directly to a file.

s = im2html(magic(10),[]);

im2html(magic(10),[],'OutputFile','magic_table.html')

Give im2html a try. Comment here (or on the File Exchange page) if you find a good use for it, or if you have ideas about making it better.

Next time I'll go into some of the details about how im2html works, including the use of raw HTML in your publishable MATLAB scripts, as well as an obscure thing in the Image Processing Toolbox called imagemodel.

Get the MATLAB code

Published with MATLAB® 7.14

Published with MATLAB® 7.14

I've written often here about various computational and visualization techniques involving labeling connected components in binary images. Sometimes I use the function label2rgb to convert a label matrix into a color image with a different color assigned to each label.

Here's an example.

bw = imread('http://blogs.mathworks.com/images/steve/2012/rice-bw.png');
imshow(bw)

Now compute the connected components and the corresponding label matrix.

cc = bwconncomp(bw);
L = labelmatrix(cc);

In the label matrix, each foreground object in the original binary image is assigned a unique positive integer. Here, for instance, is how to display the tenth object.

imshow(L == 10)

Use the function label2rgb to "colorize" the label matrix.

rgb = label2rgb(L);
imshow(rgb)

That's a nice effect, but because of the way the colors are assigned, object near each other tend to have very similar colors. It might be better sometimes to assign the colors differently. That's what the 'shuffle' argument to label2rgb is for.

Here's the full syntax including 'shuffle':

rgb = label2rgb(L,map,zerocolor,'shuffle');

zerocolor is a three-element vector specifying what color is used for the background pixels.

Let's try it with the jet colormap and a light gray background color.

rgb = label2rgb(L,'jet',[.7 .7 .7],'shuffle');
imshow(rgb)

Now it's easier to see where two objects might be actually touching and so receive the same label. Let's zoom in closer to see:

xlim([128 185]);
ylim([5 62]);

I hope you find this useful!

Get the MATLAB code

Published with MATLAB® 7.14

On my list of potential blog topics today I saw just this cryptic item labeled dftmtx. Hmm, the MATLAB dftmtx function. But have I written about this function before? I better double-check by searching the old blog postings:

Oh, I forgot about Loren's post! She showed how the discrete Fourier transform, which MATLAB users normally compute by calling fft, can also be computed via a matrix multiply.

x = rand(1000,1);
X = fft(x);

T = dftmtx(1000);
X2 = T*x;

max(abs(X2(:) - X(:)))

ans =

   3.9790e-13

The difference is just floating-point round-off error.

My old post talked about the value of having an independent computation method available when you are testing your algorithm.

So today let's do something a little different. Let's compare the performance of computing the discrete Fourier transform using the DFT matrix versus using the fast Fourier transform.

To help with the timing, I'm going to use a function I wrote called timeit that you can download from the MATLAB Central File Exchange.

clear
n = 100:50:3000;
for k = 1:length(n)
    nk = n(k);
    x = rand(nk,1);
    T = dftmtx(nk);

    f = @() T*x;
    g = @() fft(x);

    times_f(k) = timeit(f);
    times_g(k) = timeit(g);
end

plot(n,times_f,n,times_g)
legend({'Using T*x', 'Using fft(x)'})

That's a pretty dramatic difference. The blue curve, showing the computation time using T*x, is an n^2 curve. Compared to that, the green curve, showing the computation time using fft(x), is so low that you can hardly see it.

Let's expand the y axis.

ylim([0 0.0005])

The lower green curve is an n*log(n) curve. The dramatic difference between that and the n^2 curve is why everyone got so excited when the fast Fourier transform algorithm was invented (or re-invented) a few decades ago.

Get the MATLAB code

Published with MATLAB® 7.14