## Extract outline of polygons defined by gridded mask

### Kelly Kearney (view profile)

on 30 Apr 2013
Latest activity Edited by Cedric Wannaz

### Cedric Wannaz (view profile)

on 20 Sep 2017
I'm attempting to extract the polygons defined by a gridded mask, where the mask is defined pcolor-style, i.e. by the grid cell edges, and not necessarily rectilinear.
The bwboundaries function gets close to what I'm looking for, but not quite, since it does its calculations based on the center of each grid cell, rather than the edges. This causes the outline to be slightly shifted, and it leaves out one-grid-cell objects.
The following shows the discrepancy between an example mask and the outline I'm currently able to extract. Any of the image processing gurus know of a quick and easy way to extract the actual outline of all the light regions?
r = linspace(1,2,ny+1);
theta = linspace(3*pi/4,pi/4,nx+1);
[r,theta] = meshgrid(r, theta);
[xr, yr] = pol2cart(theta, r);
colormap(gray);
axis equal;
set(gca, 'clim', [-1 1]);
[xout, yout] = deal(cell(size(b1)));
for ib = 1:length(b1)
idx = sub2ind(size(xr), b1{ib}(:,1), b1{ib}(:,2));
xout{ib} = xr(idx);
yout{ib} = yr(idx);
end
hold on;
hl = cellfun(@(x,y) plot(x,y, 'b'), xout, yout);
set(hl, 'linewidth', 2);

on 30 Apr 2013

Kelly Kearney

### Kelly Kearney (view profile)

on 30 Apr 2013
Hmm, a similar issue, but my polygons are not necessarily convex, and that solution still cuts corners (in the literal sense). I was hoping to avoid the manual polygon-combo calculation, since it seems like overkill, but perhaps that's the best solution. I'll be working on that unless another answer shows me a shortcut...

### Kelly Kearney (view profile)

on 30 Apr 2013
Edited by Kelly Kearney

### Kelly Kearney (view profile)

on 30 Apr 2013

I decided to go ahead and try the not-so-quick method that first came to mind: actually calculating the polygon intersection of all connected grid cells.
% Same setup as before
r = linspace(1,2,ny+1);
theta = linspace(3*pi/4,pi/4,nx+1);
[r,theta] = meshgrid(r, theta);
[xr, yr] = pol2cart(theta, r);
% Calculate outline of each grid cell
xv{ir,ic} = [xr(ir,ic) xr(ir,ic+1) xr(ir+1,ic+1) xr(ir+1,ic) xr(ir,ic)]';
yv{ir,ic} = [yr(ir,ic) yr(ir,ic+1) yr(ir+1,ic+1) yr(ir+1,ic) yr(ir,ic)]';
end
end
% Use polybool to combine adjacent grid-cell polygons
hold on;
npoly = max(l(:));
[xpoly, ypoly] = deal(cell(npoly,1));
for ii = 1:npoly
idx = find(l == ii);
if length(idx) == 1
xpoly{ii} = xv{idx};
ypoly{ii} = yv{idx};
else
[xpoly{ii}, ypoly{ii}] = polybool('+', xv{idx(1)}, yv{idx(1)}, xv{idx(2)}, yv{idx(2)});
for ip = 3:length(idx)
[xpoly{ii}, ypoly{ii}] = polybool('+', xpoly{ii}, ypoly{ii}, xv{idx(ip)}, yv{idx(ip)});
end
end
end
This solves the problem, but as expected it's relatively slow (particularly compared to bwboundaries). And it will get even slower for more complicated grids where I may have to verify that I've defined the cell outlines in a clockwise orientation prior to calling polybool. I could probably speed things up by calling gpcmex directly (the engine under polybool), but still... polybool just seems like an overly-powerful tool to be throwing at a problem like this, and I'm sure there must be more efficient algorithms out there somewhere.

Kelly Kearney

### Kelly Kearney (view profile)

on 1 May 2013
In case anyone besides myself is interested, here's the solution I finally settled on. It basically collects the coordinates of all masked-cell edges, throws out any that are repeated more than once, and then joins the remaining edges together into polygons. The polygon-merge is still a bit slow, but it works for my purposes.
%
%
% This function finds the outline of the polygons defined by a gridded
% mask, similar to bwboundaries except that the grid is defined by its edge
% coordinates in pcolor-style. This requires that the edge definitions are
% one cell longer than the mask itself.
%
% Input variables:
%
% mx: n x 1 or m x n array defining the edges of each cell in the x
% direction
%
% my: m x 1 or m x n array defining the edges of each cell in the y
% direction
%
% mask: m-1 x n-1 array, with 1s indicating the mask region and 0 the
% background
%
% Output variables:
%
% xpoly: x coordinates of mask polygon, NaN-delimted between segments
%
% ypoly: y coordinates of mask polygon, NaN-delimted between segments
% Check input
if isvector(mx) & isvector(my)
[mx, my] = meshgrid(mx, my);
end
if ~isequal(size(mx), size(my), [nrow+1 ncol+1])
error('mx and my must be one row and column larger than mask');
end
% Edges of each cell
[xv, yv] = deal(zeros(nrow,ncol,5));
for ir = 1:nrow
for ic = 1:ncol
xv(ir,ic,:) = [mx(ir,ic) mx(ir,ic+1) mx(ir+1,ic+1) mx(ir+1,ic) mx(ir,ic)]';
yv(ir,ic,:) = [my(ir,ic) my(ir,ic+1) my(ir+1,ic+1) my(ir+1,ic) my(ir,ic)]';
end
end
xv = reshape(xv, [], 5);
yv = reshape(yv, [], 5);
xedge = xv(idx,:);
yedge = yv(idx,:);
xyedge = [...
xedge(:,1:2) yedge(:,1:2)
xedge(:,2:3) yedge(:,2:3)
xedge(:,3:4) yedge(:,3:4)
xedge(:,4:5) yedge(:,4:5)];
% Orient all edges the same way, so tracing direction doesn't
% matter
isokay = xyedge(:,1) < xyedge(:,2) | ...
(xyedge(:,1) == xyedge(:,2) & xyedge(:,3) < xyedge(:,4));
xyedge(~isokay,:) = xyedge(~isokay, [2 1 4 3]);
% Locate the edges only repeated once; these are the perimeter edges
[unqedge, nedge] = consolidator(xyedge, [], 'count');
xunq = unqedge(nedge==1, 1:2);
yunq = unqedge(nedge==1, 3:4);
xunq = [xunq'; nan(1,size(xunq,1))];
yunq = [yunq'; nan(1,size(yunq,1))];
% Merge line segments into one or more polygons
[xpoly, ypoly] = polymerge(xunq(:), yunq(:));

on 20 Sep 2017

I just ran into this same issue. My mask contains only a few grid cells and I needed an exact outline of the grid cells in my mask, preserving the 90 degree corners. Below, each black dot indicates the center of a grid cell. The true values in my mask are outlined by red circles.
In my first attempt to outline the red cells I wasn't thinking about the corners, so I tried simple contouring like this:
But as you can see, that doesn't adequately capture the corners of each grid cell. The simplest solution I've found comes via Wolfgang Schwanghart's TopoToolbox.
% Create a grid object that links x,y locations to pixel centers of the mask:
% Get the outline of the mask:
[~,xi,yi] = GRIDobj2polygon(G);
% Plot the outline in blue:
plot(xi,yi,'b')

Cedric Wannaz

### Cedric Wannaz (view profile)

on 20 Sep 2017
Another option if you have a license for ArcGIS is to build a small wrapper, and then you have access to the full toolbox. The only limitation is that the geoprocessor (or arcpy) must check the license (which may require an access to a license server). If you have many operations to perform iteratively, this means that you have to extend the mechanism so you can pass a stack of operations to your wrapper.