The developer (Brendan) who worked on this function came to see me sometime last year to find out how the 'ConvexImage' measurement offered by regionprops was computed so that he could use the same procedure for bwconvhull. I believed that I knew the answer off the top of my head, so without looking at the code I rattled off the following steps:
A few days later Brendan came back to tell me that, although my description was clear, the code that I wrote ten years ago for regionprops actually does something else.
I looked at the code and realized that Brendan was right, and I started to remember that I had actually made this same mistake many years before.
In fact, I did originally implement 'ConvexImage' in regionprops using the procedure outlined above. Before it shipped, though, I discovered (and fixed) a big problem with it.
Let me show you the problem using a small example. Here's a tiny binary image.
bw = [... 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 ]; imshow(bw, 'InitialMagnification', 'fit')
Now here's the convex image as computed by bwconvhull. The "filled-in" pixels are shown in light gray.
bwch = bwconvhull(bw) hold on h = imshow(bwch); set(h, 'AlphaData', 0.8); hold off
bwch = 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 1 1 0 0 0 0 0 1 1 1 1 1 0 0 0 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0
Now I'll graphically illustrate the computational procedure above. First, compute the set of corner locations for all the foreground pixels.
[y,x] = find(bw); dx = [-0.5 -0.5 0.5 0.5]; dy = [-0.5 0.5 -0.5 0.5]; x_corners = bsxfun(@plus, x, dx); y_corners = bsxfun(@plus, y, dy); x_corners = x_corners(:); y_corners = y_corners(:);
Visualize step 1 by superimposing those corner locations on the input image.
imshow(bw, 'InitialMagnification', 'fit') hold on plot(x_corners, y_corners, 'og') hold off
Step 2: Compute the convex hull of all the corner points. Visualize by superimposing the convex hull in red.
k = convhull(x_corners, y_corners); x_hull = x_corners(k); y_hull = y_corners(k); hold on plot(x_hull, y_hull, 'r', 'LineWidth', 4) hold off
Step 3: Use poly2mask to convert the convex hull to a binary image mask.
mask = poly2mask(x_hull, y_hull, 6, 9); imshow(bw, 'InitialMagnification', 'fit') hold on h = imshow(mask); set(h, 'AlphaData', 0.8) plot(x_hull, y_hull, 'r', 'LineWidth', 4) hold off
You can see that there are extra pixels along the diagonal edges that got put into the mask. That's not good. If you repeat the operation, those diagonals will keep filling out until you're left with a rectangle. That's not supposed to happen with repeated application of the convex hull.
The reason this is happening is that the convex hull goes exactly through the center of pixels that are along the diagonal but "outside" the original set of foreground pixels. Because of the tie-breaking rules applied by poly2mask, all (in this case) of those outside pixels got included in the output mask.
The solution that I settled on ten years ago, and which is now also used in bwconvhull, was to modify the first step of the procedure. Instead of collecting the set of corner points of each foreground pixel, the correct procedure collects points along the center of each edge of each foreground pixel.
Here's what that looks like.
dx = [0.0 0.0 0.5 -0.5]; dy = [0.5 -0.5 0.0 0.0]; x_edges = bsxfun(@plus, x, dx); y_edges = bsxfun(@plus, y, dy); x_edges = x_edges(:); y_edges = y_edges(:); k = convhull(x_edges, y_edges); x_hull = x_edges(k); y_hull = y_edges(k); imshow(bw, 'InitialMagnification', 'fit') hold on plot(x_edges, y_edges, 'og') plot(x_hull, y_hull, 'r', 'LineWidth', 4) hold off
mask = poly2mask(x_hull, y_hull, 6, 9); imshow(mask, 'InitialMagnification', 'fit') hold on plot(x_hull, y_hull, 'r', 'LineWidth', 4) plot(x_edges, y_edges, 'og') hold off
I have been fooled more often than I would like to admit by sometimes nonintuitive digital image geometry. For you image processing algorithm people out there, I hope seeing these pictures and techniques will help you avoid similar pitfalls someday.
Get the MATLAB code
Published with MATLAB® 7.13
Comments are closed.
9 CommentsOldest to Newest
interesting. I’m keen to know if there’s a proof that this does what you expect: that is, whether bwconvhull(bwconvhull(X)) is equal to bwconvhull(X)
Andrew—I have not attempted such a proof.
Very neat and obvious function to have. It still only appears to accept 2d images, though. Here’s to hoping for a 3d or n-d implementation for it in the future so we’re not stuck with vert2con!
(Also would be nice in regionprops!)
Sean—Thanks for the suggestion.
Sean—Can you provide us with one or examples of specific applications or use cases for your request?
Sure. For my MS research one small task was to identify glass beads in CT scans. The segmentation based on standard deviation worked well, but failed in the presence of ring/beam artifacts. Since we knew the beads were ideally convex (occasionally there would be two in a figure eight shape but those could be ignored), a simple way to fill the ring artifact troughs cut through the beads was to calculate each bead’s convex hull and automatically fill it, basically using what appears to be the ‘objects’ option in bwconvhull.
%Connect components, part 1, bounding box/image to fill with convex hull CC = bwconncomp(Xbeads); RP = regionprops(CC,'Image','boundingbox'); RP = fillConvexHull(RP); %Replace Beads in Image with Convexhull for ii = 1:length(RP) BBsz = ceil(RP(ii).BoundingBox); Xbeads(BBsz(2):BBsz(2)+BBsz(5)-1,BBsz(1):BBsz(1)+BBsz(4)-1,BBsz(3):BBsz(3)+BBsz(6)-1) = RP(ii).Image|Xbeads(BBsz(2):BBsz(2)+BBsz(5)-1,BBsz(1):BBsz(1)+BBsz(4)-1,BBsz(3):BBsz(3)+BBsz(6)-1); end
with fillConvexHull() being
function RP = fillConvexHull(RP) %Fill the convex hull of each bead. for ii = 1:length(RP) idx = find(RP(ii).Image); sz = size(RP(ii).Image); [r c p] = ind2sub(sz,idx); [A,b] = vert2con([r c p]); [rr cc pp]=ndgrid(1:sz(1),1:sz(2),1:sz(3)); V=[rr(:) cc(:) pp(:)]'; V=all(bsxfun(@le,A*V,b)); RP(ii).Image=reshape(V,sz); end end
Sean—Thanks, that’s helpful.
I was playing around with convex hull algorithms last month. Your idea here is very good, I hope you don’t mind me taking it!
But, and answering Andrew’s question, there seems to be no way of computing a convex hull in a discrete grid such that the function becomes idempotent.
bw = [... 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 1 0 0 1 1 0 0 0 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 ]; bw = bwconvhull(bw); isequal(bw,bwconvhull(bw))
I’ve tried with using only the centre of each pixel, and even then some straight lines keep moving outward when repeatedly taking the convex hull. It must be a property of discrete lines, and I bet someone’s already written a proof… :)
Cris—Nice to hear from you. I’m glad you like the idea, and I would be delighted if you made good use of it.
From my experience with geometry algorithms in image processing, it doesn’t surprise me that you quickly found an example proving that bwconvhull is not idempotent. I guess it would challenging (and perhaps impossible, as you suggest) to do it.
My feelings about geometrical image processing algorithms are heavily influenced by reading dozens of papers on image thinning in years past. As far as I can tell, all papers about image thinning have the following form: “We demonstrate that these other N thinning methods to be wrong in several cases, and we present a new method that is correct.”