Smallest mask enclosing a polygon
10 次查看(过去 30 天)
显示 更早的评论
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 :
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1467211/image.png)
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 个评论
采纳的回答
Bruno Luong
2023-8-30
编辑:Bruno Luong
2023-8-30
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.
3 个评论
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.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1468131/image.png)
Bruno Luong
2023-8-30
编辑:Bruno Luong
2023-8-30
Somewhat more robust version, since interp1 can throw an error in some rare circumstances.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1468496/image.png)
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
更多回答(3 个)
Image Analyst
2023-8-30
Perhaps you'd be interested in the attached paper on "Minimum Perimeter Polygon"
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1468431/image.png)
0 个评论
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;
% mask = zeros(size(X));
[in,on] = inpolygon(X,Y,shape(1,:),shape(2,:)) ;
mask(i,j)=in|on;
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 个评论
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 = load('example.mat')
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)
binaryImage = false(rows, columns);
[rowss, columnss] = size(xs)
% xs and ys are broken into groups separated by NaN. So find each groups.
nanLocations = [1; find(isnan(xs)); rowss]
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
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!