Asked by Kelly Kearney
on 30 Apr 2013

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 by Sean de Wolski
on 30 Apr 2013

This reminded me of this blog post by Steve:

Kelly Kearney
on 30 Apr 2013

Sign in to comment.

Answer by Kelly Kearney
on 30 Apr 2013

Edited by Kelly Kearney
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

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

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

Sign in to comment.

Answer by Chad Greene
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:

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 Wannaz
on 20 Sep 2017

Sign in to comment.

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.