# Smallest mask enclosing a polygon

6 ビュー (過去 30 日間)
Philippe EAR 2023 年 8 月 29 日

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
##### 2 件のコメントなしを表示なしを非表示
Dyuman Joshi 2023 年 8 月 29 日

"Is there a function for that "
Philippe EAR 2023 年 8 月 29 日

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 2023 年 8 月 30 日

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 件のコメント1 件の古いコメントを表示1 件の古いコメントを非表示
Bruno Luong 2023 年 8 月 30 日
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 2023 年 8 月 30 日

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

サインインしてコメントする。

### その他の回答 (3 件)

Image Analyst 2023 年 8 月 30 日
Perhaps you'd be interested in the attached paper on "Minimum Perimeter Polygon" ##### 0 件のコメント-2 件の古いコメントを表示-2 件の古いコメントを非表示

サインインしてコメントする。

KSSV 2023 年 8 月 29 日
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,:)) ;
##### 1 件のコメント-1 件の古いコメントを表示-1 件の古いコメントを非表示
Philippe EAR 2023 年 8 月 29 日
Thanks for the tip on the loop.
I cannot increase the resolution, my longitude and latitude meshgrid are fixed.

サインインしてコメントする。

Image Analyst 2023 年 8 月 29 日
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
##### 5 件のコメント3 件の古いコメントを表示3 件の古いコメントを非表示
Philippe EAR 2023 年 8 月 29 日
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 2023 年 8 月 29 日
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 サインインしてコメントする。

### カテゴリ

Help Center および File ExchangeColormaps についてさらに検索

R2022b

### Community Treasure Hunt

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

Start Hunting!