Contents

I’ve been thinking a lot about how to generate and efficiently store paths–and then reconstruct binary images from them. Path simplification is important in many fields, and after some Googling, I found my way to the “Douglas-Peucker” algorithm. (It goes by other names as well, but that seems to be the most common.)

It turns out that Douglas-Peucker is used extensively in cartography to simplify the representation of geographical borders. Armed with that information, I came back to MATLAB’s doc-search engine and, by typing “Douglas Peucker,” I found my way to the Mapping Toolbox implementation of reducem. From the doc:

The reducem function implements a powerful line simplification algorithm (known as Douglas-Peucker) that intelligently selects and deletes visually redundant points.

For those of you (us) who have the Mapping Toolbox, you’re good to go. Despite the fact that it was implemented for a non-Cartesian geometry, it works nicely, as I will demonstrate. For others, a quick check of the File Exchange reveals four implementations of this particular algorithm:

Let’s create a sample path, and consider each of them (along with reducem) in turn:

Sample Path

pgon = nsidedpoly(5,'Center',[150 150],'SideLength',100);
verts = pgon.Vertices;
verts = [verts;verts(1,:)];
verts = [verts(1:2:end,:);verts(2:2:end,:)];
mask = poly2mask(verts(:,2),verts(:,1),300,300);
mask = imclose(mask,strel('disk',4));
mask = imfill(mask,'holes');
figure
subplot(1,2,1)
imshow(mask)
title('BW Star')
perim = bwboundaries(mask);
perim = perim{1};
subplot(1,2,2)
plot(perim(:,2),perim(:,1),'b-')
title('Extracted Boundary')
axis equal
set(gca,'xlim',[0,300]+0.5, 'ylim',[0,300]+0.5)

If you run that code, you’ll see that variable perim has 547 [row, col] coordinate pairs. Most of them are redundant for our purposes, and we would like to be able to discard them.

Enter the Douglas-Peucker algorithm:

x = perim(:,2);
y = perim(:,1);
simplified = DouglasPeucker([x';y'],0)

Wait…what’s going on there? Turns out that Reza’s implementation considers the distance between the first point and the last. That isn’t very useful for a closed curve! I’m going to move on to…

[list_pts,pnts]= douglas_peucker([x,y], 0);
size(pnts)
ans =
     2     2

Here again, the results match Reza’s (the two returned points are the point on the “upper left” prong, and are coincident). Ranadeer’s implementation also relies on distances from the endpoints of the curve. (Note that changing the tolerance in either Reza’s or Ranadeer’s implementation has no effect.) Let’s try…

simplified = DouglasPeucker([x,y],10);

Beautiful! That’s what I’m talking about! Conveniently, Ligong coded to allow an additional input flag to show the data:

Nice!

Next up, Line Simplification, by Wolfgang Schwanghart:

[ps,ix] = dpsimplify([x,y],10);

Again, that’s lovely! (Note that Wolfgang also outputs idx–an index of the points that are retained. That could be useful!)

Finally, consider the Mapping Toolbox implementation. Even though it is implemented for geospatial (latitude/longitude) data, it works beautifully!

[xReduced,yReduced] = reducem(perim(:,2),perim(:,1),10);
plot(xReduced,yReduced,'r.','Markersize',14)

Mask reconstruction

So clearly we have some great options, even for simplifying closed curves. Once you’ve got that reduced set of vertices, note that you can readily reconstruct the polygon using poly2mask():

mask2 = poly2mask(ps(:,1),ps(:,2), 300, 300);

It is important to note that, while those binary images are nearly identical, they are not exactly identical. Why? Steve (in one of his stellar blog posts on image processing) discusses that issue. It turns out that we can get identical results, if we upsample the path, with nearest-neighbor interpolation, by a factor of 3!

As always, I welcome your thoughts and comments.

Get the MATLAB code

Published with MATLAB® R2018a