Brett recently asked me about an image processing problem. One step in the problem involved finding out which pixels are within a certain distance of the foreground in a binary image. We experimented with a few methods, and I want to share the results with you. In particular, I want to show you a useful technique called isotropic dilation that can be achieved using the distance transform.
Let's start with a small 400-by-400 test image that has only a single foreground pixel in the center.
bw = false(400, 400); bw(200, 200) = 1; imshow(bw)
Suppose we want to know which pixels are within 25 units of the foreground pixel? One common technique is to perform dilation with a disk.
se = strel('disk', 25);
bw2 = imdilate(bw, se);This has a problem, though. If you read the documentation for strel carefully, you'll see that by default it computes a structuring element that is an approximation of a disk. You can see the effect of the approximation by looking at the output of the dilation:
imshow(bw2)
You can avoid the approximation by passing in 0 as a third input argument to strel:
se2 = strel('disk', 25, 0);
bw3 = imdilate(bw, se2);
imshow(bw3)
So why does strel use a disk approximation by default? Because it's usually a lot faster. Let's check it out on a larger test image (750-by-1500) with a lot of pixels in the foreground.
url = 'http://blogs.mathworks.com/images/steve/2010/binary_benchmark.png'; bw2 = imread(url); imshow(bw2, 'InitialMagnification', 25)
f = @() imdilate(bw2, se) time_with_approximate_disk = timeit(f)
f =
@()imdilate(bw2,se)
time_with_approximate_disk =
0.0292
g = @() imdilate(bw2, se2) time_with_exact_disk = timeit(g)
g =
@()imdilate(bw2,se2)
time_with_exact_disk =
1.4536
(You can download timeit from the MATLAB Central File Exchange.)
An alternative to using imdilate is to use bwdist to compute the distance transform and then threshold the result. The distance transform computes, for each pixel, the distance between that pixel and the nearest foreground pixel. Here's how it would work:
bw4 = bwdist(bw) <= 25; imshow(bw4)
The function bwdist uses a fast algorithm, so this works a lot faster than imdilate with an exact disk.
h = @() bwdist(bw2) <= 25; time_using_distance_transform = timeit(h)
time_using_distance_transform =
0.0687
This technique of thresholding the distance transform is sometimes called isotropic dilation.
Have you made creative use of the distance transform in your work? Then post a comment here.
Get
the MATLAB code
Published with MATLAB® 7.10
Not sure if this really is creative. We have used it to find the pixel of feature-B closest to feature-A of the image. We(my biologist friend and me) were studying beating pattern of the tail of sea-urchin sperm. After segmenting the head and the tail, we needed to locate the point on the tail that was closest to head. I computed distance transform with head as the foreground pixel, masked it with the segmentation of the tail, and sorted the masked distance transform to find the pixel!
I use bwdist quite a bit with medical images. Things like finding bits of the body that can be expected to be near the skin/lungs etc.
One question I have is the following:
bwdist uses (as you state) isotropic dilation. This is fine for 2D medical images since most images are in fact isotropic. If I want to convert bwdist output to true spatial distance, it’s a simple case of scaling by the pixel-to-mm ratio.
When considering 3D volumes however, you’ll often get an interslice spacing that’s different to the slice pixel size. If I want to get a true spatial bwdist, I need to somehow create a non-isotropic dilation.
It would be great to be able to supply pixel-spacing as extra variables to bwdist (see the gradient() function) to facilitate this.
Im am working on a project for detecting several features of a training exercice in realtime. The code is implemented in C#, but the image processing toolbox I’m using doesn’t support any skeletonization. Therefore I am using the distance transform in combination with an LoG filter to “approximate” the skeletonization. The results are quite good.
I’ve used bwdist to speed up any erosion/dilation step that requires manual interaction to set. For example, you can calculate bwdist once, then use a sliderbar to set the threshold/distance very quickly.
Jeff—Nice! I hadn’t thought of that.
Sven—There might be a way to accomplish this by using a bwdist helper function (private/ddist) directly. I’ll think about it and maybe do a blog post on it.
Thanks Steve – at the moment I will usually use a MEX nearest-neighbour finder from the file exchange and do the following:
1. Collect the XYZ locations of all ON and OFF b/w voxels.
2. Find the nearest ON neighbour index of the OFF list.
3. Calculate the distance from each OFF voxel to its nearest ON pair
As you can imagine, this gets pretty memory-intensive so I need to convert to singles wherever possible, and wrap step 3 in a try/catch statement that tries to calculate all distances in one command, but resorts to looping over each pixel if a memory error occurs.
As you suggest, it looks like bwdist internally stores some weights before passing to ddist. Thanks for the reply, this looks promising to me :)
Sven—Yes, you caught my hint correctly. :-) Just be aware that the two-pass algorithm used by private/ddist with the weights vector isn’t an exact algorithm.
Steve, Thanks for blogging about this—great post, great comments!
I found it the bwdistance transform approach exceptionally useful, and am using it precisely as Jeff described. In your introductory sentence, you said that I was looking for pixels within a certain distance of a foreground binary object. More precisely, I was looking for _regions_ within a certain distance of other regions; for the problem at hand, I wanted to treat “close” ROIs as non-unique. In other words, in a blob analysis problem, any two ROIs that are discontiguous (using whatever measure of connectivity) are typical treated as separate blobs. I wanted to be able to detect cell colonies in a dish, and to say that any colony within _minSep_ pixels of another should be considered (by |regionprops|) to be a non-distinct from its neighbor.
The first problem was determining _minSep_ in this regard; a slider UIcontrol made short work of that. Behind the callback of that slider, I converted the output of BWDIST to binary by comparing it to minSep/2 (read from the slider), and then BWLABELed that binary matrix. The result properly labels the regions, treating as distinct any ROIs that are farther apart than minSep, and as non-distinct any that are within that distance. However, that dilated the original ROIs, so I followed the BWLABEL step with a logical mask, to indicate that any pixel that was black in the original binary image should be black in the new binary image as well:
D = bwdist(bw); bw2 = (D <= minSep/2); L = bwlabel(bw2); L(~bw) = 0;The result is a “minimum separation slider,” quite useful for the cell colony problem! I had originally implemented this with a standard morphological dilation, but this was approach was _much_ faster. So thanks for this!
How ya doing there Steve
I have a binary image with lots of columns of objects in it, I am trying to find distances between the various columns of objects closest to each other, from the center points of each. Is there a Matlab function to do this ?????? I have already calculated the center points of each column of object in a seperate image.
I noticed in some of your above examples you try similar stuff.
Regards Jack