Extract outline of polygons defined by gridded mask

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?

 I = imread('rice.png');
 mask = im2bw(I, graythresh(I));
 mask = mask(1:100,1:100);
 [ny,nx] = size(mask);

 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);

 hp = pcolor(xr,yr,[mask nan(ny,1); nan(1,nx+1)]);
 colormap(gray);
 axis equal;
 set(gca, 'clim', [-1 1]);

 b1 = bwboundaries(mask, 4);

 [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);

Answer 1

[Sean de Wolski]

This reminded me of this blog post by Steve:

http://blogs.mathworks.com/steve/2011/10/04/binary-image-convex-hull-algorithm-notes/

[Kelly Kearney]

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...

Answer 2

[Kelly Kearney]

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

 I = imread('rice.png');
 mask = im2bw(I, graythresh(I));
 mask = mask(1:100,1:100);
 [ny,nx] = size(mask);

 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

 for ir = 1:size(mask,1)
     for ic = 1:size(mask,2)
         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

 l = bwlabel(mask, 4);
 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]

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.

function [xpoly,ypoly] = mask2poly(mx, my, mask)
%MASK2POLY Calculates polygon outline of a gridded mask
%
% [x,y] = mask2poly(mx, my, mask)
%
% 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

% Copyright 2013 Kelly Kearney

% Check input

if isvector(mx) & isvector(my)
    [mx, my] = meshgrid(mx, my);
end

[nrow, ncol] = size(mask);

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);

idx = find(mask);
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(:));

Answer 3

[Chad Greene]

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:

contour(x,y,double(mask),[0.5 0.5],'r')

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: 
G = GRIDobj(x,y,mask);

% Get the outline of the mask: 
[~,xi,yi] = GRIDobj2polygon(G);

% Plot the outline in blue: 
plot(xi,yi,'b')

[Cedric]

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.