# Smallest mask enclosing a polygon

9 views (last 30 days)
Philippe EAR on 29 Aug 2023
Edited: Bruno Luong on 30 Aug 2023
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;
for i =1:size(X,1)
for j = 1:size(X,2)
if inpolygon(X(i,j),Y(i,j),shape(1,:),shape(2,:))
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 on 29 Aug 2023
Moved: Dyuman Joshi on 30 Aug 2023
"Is there a function for that "
Philippe EAR on 29 Aug 2023
Moved: Dyuman Joshi on 30 Aug 2023
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

Bruno Luong on 30 Aug 2023
Edited: Bruno Luong on 30 Aug 2023
To me there is no other way than detecting the border cross the edges of each "pixel"
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));
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);
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)
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)
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.
##### 3 CommentsShow 1 older commentHide 1 older comment
Bruno Luong on 30 Aug 2023
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 on 30 Aug 2023
Edited: Bruno Luong on 30 Aug 2023
Somewhat more robust version, since interp1 can throw an error in some rare circumstances.
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;
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);
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)
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)
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

Image Analyst on 30 Aug 2023
Perhaps you'd be interested in the attached paper on "Minimum Perimeter Polygon"

KSSV on 29 Aug 2023
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;
[in,on] = inpolygon(X,Y,shape(1,:),shape(2,:)) ;
Philippe EAR on 29 Aug 2023
Thanks for the tip on the loop.
I cannot increase the resolution, my longitude and latitude meshgrid are fixed.

Image Analyst on 29 Aug 2023
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 on 29 Aug 2023
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));
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 on 29 Aug 2023
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 = 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

### Categories

Find more on Image Arithmetic in Help Center and File Exchange

R2022b

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!