Smallest mask enclosing a polygon

I want to find the smallest mask that contains a given polygon (for example the shape of a country).

X = longitude;
Y = latitude;
[Y,X] = meshgrid(Y,X);
shape = some_polygon;
mask = zeros(size(X));
for i =1:size(X,1)
    for j = 1:size(X,2)
        if inpolygon(X(i,j),Y(i,j),shape(1,:),shape(2,:))
            mask(i,j)=1;
        end
    end
end

X and Y are lon/lat grids.

This would only give me the coordinates that are in or on the polygon, so if I were to plot the given mask and the polygon overlapped, I may have areas inside the polygon that are not filled.

Example of what that would give me :

I'd like to have the (smallest) mask that would completly enclose the polygon and thus not have any white space.

Is there a function for that

[Dyuman Joshi]

"Is there a function for that "

Try convhull.

[Philippe EAR]

I would need to be able to put constraints of which points convhull can use to create the englobing "subgrid", and it may not be convex either

Answer 1

[Bruno Luong]

To me there is no other way than detecting the border cross the edges of each "pixel"

load('example.mat')

P = example.shape;
X = example.X;
Y = example.Y;
x = X(:,1);
y = Y(1,:);

dx = mean(diff(x));
dy = mean(diff(y));

xe = [x(1)-dx/2; (x(1:end-1)+x(2:end))/2; x(end)+dx/2];
ye = [y(1)-dy/2, (y(1:end-1)+y(2:end))/2, y(end)+dy/2];
nx = length(xe);
ny = length(ye);

Pi = [interp1(xe(:), (1:nx)', P(:,1)), ...
      interp1(ye(:), (1:ny)', P(:,2))];

n = size(Pi,1);
k = 1;
A0 = Pi(1,:);
A = A0;
in = 2*inpolygon(X,Y,P(:,1),P(:,2));
bdrmask = 1;
while k < n
    k = k + 1;
    B = Pi(k,:);
    newcomponent = any(isnan(B));
    if newcomponent
        % wrap around
        B = A0;
    end
    xA = A(1); xB = B(1);
    yA = A(2); yB = B(2);
    in(floor(xA),floor(yA)) = bdrmask;
    in(floor(xB),floor(yB)) = bdrmask;
    if xB < xA
        xi = ceil(xB):floor(xA);
    else
        xi = ceil(xA):floor(xB);
    end
    if ~isempty(xi)
        yi = interp1([xA xB], [yA yB], xi);
        for i=1:length(xi)
            in(xi(i)+[-1 0],floor(yi(i))) = bdrmask;
        end
    end

    if yB < yA
        yi = ceil(yB):floor(yA);
    else
        yi = ceil(yA):floor(yB);
    end
    if ~isempty(yi)
        xi = interp1([yA yB], [xA xB], yi);
        for i=1:length(yi)
            in(floor(xi(i)),yi(i)+[-1 0]) = bdrmask;
        end
    end

    if newcomponent
        k = k+1;
        if k > n
            break
        end
        B = Pi(k,:);
        A0 = B;
    end
    A = B;
end

figure('position',[258 209 582 524]);
imagesc(x,y,double(in'))
set(gca,'Ydir','normal')
hold on
plot(P(:,1),P(:,2),'k','linewidth',0.5)
yline = ye([1 end]);
for i=1:nx
    plot(xe(i)+[0 0], yline, 'color', 0.5+[0 0 0]);
end
xline = xe([1 end]);
for i=1:ny
    plot(xline, ye(i)+[0 0], 'color', 0.5+[0 0 0]);
end

colormap jet

WARNING: the code supposes the grid contain the region(s) with at least one pixel margin.

It is also assumes te grid is more or less regular since the polygonal border "shape" is assumed to be a line in the index space, which is NOT exactly in the coordinate space in case the coordinate grid is not uniform.

[Philippe EAR]

Thank you for the solution, it seems to work fine for my application since I'm only working on regular grids.

What do you mean by a "one pixel margin" ?

One solution I think may work if the shape is fine enough is for each vertex, find the corresponding grid cell that contains that vertex. It should give the same green shape but some edge cases of border crossing can be found, but the "forgotten" areas should be limited by how fine the shape is.

[Bruno Luong]

At a second though, no a pixel margin is not required. I was afraid of indexing -1 by

in(xi(i)+[-1 0],floor(yi(i)))

but it would never happen as xe/ye extends half pixel. So shape must be within the region (not goes beyond that).

I don't think reduce the finess of boundary can help. assuming the shape has a line as the red segment, the blue pixel must be counted but it does not contain any vertex of shape. and it can cross like that regardless the resolution of shape vs pixel.

So the only way to deal with such eventual issue is detecting pixel edge crossing for all lines of shape.

[Bruno Luong]

Somewhat more robust version, since interp1 can throw an error in some rare circumstances.

load('example.mat')

bdrmask = 1;
insidemask = 2;

P = example.shape;
X = example.X;
Y = example.Y;
x = X(:,1);
y = Y(1,:);

dx = mean(diff(x));
dy = mean(diff(y));

xe = [x(1)-dx/2; (x(1:end-1)+x(2:end))/2; x(end)+dx/2];
ye = [y(1)-dy/2, (y(1:end-1)+y(2:end))/2, y(end)+dy/2];
nx = length(xe);
ny = length(ye);

Pi = [interp1(xe(:), (1:nx)', P(:,1), 'linear'), ...
      interp1(ye(:), (1:ny)', P(:,2), 'linear')];

n = size(Pi,1);
k = 1;
A0 = Pi(1,:);
A = A0;
in = insidemask*inpolygon(X,Y,P(:,1),P(:,2));
while k < n
    k = k + 1;
    B = Pi(k,:);
    newcomponent = any(isnan(B));
    if newcomponent
        % wrap around
        B = A0;
    end
    xA = A(1); xB = B(1);
    yA = A(2); yB = B(2);
    in(floor(xA),floor(yA)) = bdrmask;
    in(floor(xB),floor(yB)) = bdrmask;
    if xB < xA
        xi = ceil(xB):floor(xA);
    else
        xi = ceil(xA):floor(xB);
    end
    if ~isempty(xi) && (xA~=xB)
        yi = floor(yA + (yB-yA)/(xB-xA)*(xi-xA));
        for i=1:length(xi)
            in(xi(i)+[-1 0], yi(i)) = bdrmask;
        end
    end

    if yB < yA
        yi = ceil(yB):floor(yA);
    else
        yi = ceil(yA):floor(yB);
    end
    if ~isempty(yi) && (yA~=yB)
        xi = floor(xA + (xB-xA)/(yB-yA)*(yi-yA));
        for i=1:length(yi)
            in(xi(i),yi(i)+[-1 0]) = bdrmask;
        end
    end

    if newcomponent
        k = k+1;
        if k > n
            break
        end
        B = Pi(k,:);
        A0 = B;
    end
    A = B;
end

close all
figure('position',[258 209 582 524]);
imagesc(x,y, in.')
set(gca,'Ydir','normal')
hold on
yline = ye([1 end]);
for i=1:nx
    plot(xe(i)+[0 0], yline, 'color', 0.5+[0 0 0]);
end
xline = xe([1 end]);
for i=1:ny
    plot(xline, ye(i)+[0 0], 'color', 0.5+[0 0 0]);
end
plot(P(:,1),P(:,2),'k','linewidth',0.5)
colormap jet

Answer 2

[Image Analyst]

Perhaps you'd be interested in the attached paper on "Minimum Perimeter Polygon"

Answer 3

[KSSV]

You need not to use a loop. If you want fine mask, increase the resolution.

X = longitude;
Y = latitude;
[Y,X] = meshgrid(Y,X);
shape = some_polygon;
% mask = zeros(size(X));
[in,on] = inpolygon(X,Y,shape(1,:),shape(2,:)) ; 
mask(i,j)=in|on;

[Philippe EAR]

Thanks for the tip on the loop.

I cannot increase the resolution, my longitude and latitude meshgrid are fixed.

Answer 4

[Image Analyst]

Not sure what you're calling mask and what you're calling polygon, but the smallest shape that would enclose the blue shape without any white space would be the blue shape itself.

If you want to smooth boundaries, there is a variety of ways to do that, like with activecontour or by simply blurring the shape and thresholding, or by using imclose

[Philippe EAR]

The polygon here would be the outline of France. I have a fixed grid of longitude/latitude of which I cannot change the resolution of.

What I would like is the smallest shape made of square (corresponding to my lon/lat resolution) that completely include the outline of France.

[Image Analyst]

OK, so the polygon is the perimeter of the blue shape, not the black line.

So the smallest shape is the blue shape itself. That is the shape that will have the smallest area. To have the shape with the smallest perimeter, use convhull.

Is your data a matrix? And do you want the perimeter of the blue shape? If so you can treat the blue France as an image and use bwboundaries to get a list of (x,y) coordinates of the perimeter of the blue, or use bwperim to get the perimeter as an image where it's 1 along the perimeter and zero everywhere else, including inside.

Attach your data in a .mat file if you can.

[Dyuman Joshi]

@Philippe EAR

You said this in the problem statement -

"I'd like to have the (smallest) mask that would completly enclose the polygon and thus not have any white space."

And you said this above -

"What I would like is the smallest shape made of square (corresponding to my lon/lat resolution) that completely include the outline of France."

We can't help you if you keep changing what you want to obtain.

What is it that you want? Clarify it properly once.

Also, it will be better to show a rough image of what you want to achieve.

[Philippe EAR]

What I call mask is a binary matrix that is True if the corresponding (x,y) coordinates (X(i,j),Y(i,j)) is needed to completely include the black outline.

The current inpolygon function would create something like this :

Where only the cells that are completely inside the black outline are selected.

But I want to find the cells that would create the smallest blue shape that includes the black outline, so something like :

Here's the code to visualize my issue with the associated .mat file:

[in,on] = inpolygon(example.X,example.Y,example.shape(:,1),example.shape(:,2));
mask = double(in|on);
surf(example.X,example.Y,mask);view(2);shading flat; hold on;
plot3(example.shape(:,1),example.shape(:,2),zeros(size(example.shape(:,1),1),1)+1000000,"-k","linewidth",2)

Here is the result, I want the yellow shape to be the smallest possible, while completely containing the black outline.

[Image Analyst]

There are 3 shapes separated by a NaN in the list.

This is what I have so far but I don't think we have enough resolution on the grid.

s = load('example.mat')

s = struct with fields:
    example: [1×1 struct]

x = s.example.X;
y = s.example.Y;
shape = s.example.shape;
xs = shape(:, 1);
ys = shape(:, 2);
subplot(1, 2, 1);
plot(xs, ys, 'r-', 'LineWidth', 2)
grid on;
axis equal;
legend
% Shift so that min ia 1 so we can make an image.
xs = xs - min(xs) + 1;
ys = ys - min(ys) + 1;
% Get a digital image
[rows, columns] = size(x)

rows = 69

columns = 49

binaryImage = false(rows, columns);
[rowss, columnss] = size(xs)

rowss = 19499

columnss = 1

% xs and ys are broken into groups separated by NaN.  So find each groups.
nanLocations = [1; find(isnan(xs)); rowss]

nanLocations = 4×1

1 16978 17000 19499

for k = 1 : length(nanLocations) - 1
    index1 = nanLocations(k) + 1;
    index2 = nanLocations(k+1) - 1;
    fprintf('Getting boundary #%d between %d and %d\n', k, index1, index2);
    thisx = xs(index1 : index2);
    thisy = ys(index1 : index2);
    binaryImage = binaryImage | poly2mask(thisx, thisy, rows, columns);
    subplot(1, 2, 2);
    imshow(binaryImage)
    axis('on', 'xy')
    drawnow
end

Getting boundary #1 between 2 and 16977
Getting boundary #2 between 16979 and 16999
Getting boundary #3 between 17001 and 19498