Help finding datapoints in x,y falling outside the boundaries of a polygon. Incongruencies between inpolygon and inpoly2, time issues, and possible solutions

5 次查看（过去 30 天）
Simone A.
Simone A. 2023-10-2
移动：Bruno Luong 2023-10-5
Dear all,
To provide a visual explanation (see below), I have a scatterplot x,y. Overlayed, I have a polygon (let's call it P). I 'simply' need to find all the datapoints (x,y) falling outside the boundaries of P.
I can achieve this by using inpolygon. This however slows down my code by ~10 seconds (please note I need to repeat the inpolygon step several times within my code). To solve for the time issue, I tried inpoly2 ( https://it.mathworks.com/matlabcentral/fileexchange/10391-inpoly-a-fast-points-in-polygon-test?s_tid=prof_contriblnk ), which performs exactly as inpolygon in most occasions. This reduces the run time by ~80% which is great!
Yet, in some circumstances inpoly2 fails to correctly identify the points enclosed within the polygon, which means that I cannot use it unless I am able to find a solution.
Please see below (find attached the variables used in this code for reproducibility):
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Method 1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
is_outside1 = ~inpolygon(xy(:, 1), xy(:, 2), p1.Vertices(:, 1), p1.Vertices(:, 2));
Idx_Out_1 = find(is_outside1);
is_outside2 = ~inpolygon(xy(:, 1), xy(:, 2), p2.Vertices(:, 1), p2.Vertices(:, 2));
Idx_Out_2 = find(is_outside2);
is_outside3 = ~inpolygon(xy(:, 1), xy(:, 2), p3.Vertices(:, 1), p3.Vertices(:, 2));
Idx_Out_3 = find(is_outside3);
is_outside4 = ~inpolygon(xy(:, 1), xy(:, 2), p4.Vertices(:, 1), p4.Vertices(:, 2));
Idx_Out_4 = find(is_outside4);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Method 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[is_outside11, ~] = inpoly2(xy, p1.Vertices);
Idx_Out_11 = find(is_outside11==0);
[is_outside22, ~] = inpoly2(xy, p2.Vertices);
Idx_Out_22 = find(is_outside22==0);
[is_outside33, ~] = inpoly2(xy, p3.Vertices);
Idx_Out_33 = find(is_outside33==0);
[is_outside44, ~] = inpoly2(xy, p4.Vertices);
Idx_Out_44 = find(is_outside44==0);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Plot %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure
subplot(2,2,1)
plot(xy(:, 1),xy(:, 2),'k.')
hold on
plot(p1,'FaceColor','none','EdgeColor','r', 'LineWidth',2)
plot(xy(Idx_Out_1, 1),xy(Idx_Out_1, 2),'xr')
plot(xy(Idx_Out_11, 1),xy(Idx_Out_11, 2),'ob')
legend('data','p1 Poly','Outside inpolygon','Outside inpoly2',Location='northwest')
title('P1')
subplot(2,2,2)
plot(xy(:, 1),xy(:, 2),'k.')
hold on
plot(p2,'FaceColor','none','EdgeColor','y', 'LineWidth',2)
plot(xy(Idx_Out_2, 1),xy(Idx_Out_2, 2),'xr')
plot(xy(Idx_Out_22, 1),xy(Idx_Out_22, 2),'ob')
legend('data','p2 Poly','Outside inpolygon','Outside inpoly2',Location='northwest')
title('P2')
subplot(2,2,3)
plot(xy(:, 1),xy(:, 2),'k.')
hold on
plot(p3,'FaceColor','none','EdgeColor','g', 'LineWidth',2)
plot(xy(Idx_Out_3, 1),xy(Idx_Out_3, 2),'xr')
plot(xy(Idx_Out_33, 1),xy(Idx_Out_33, 2),'ob')
legend('data','p3 Poly','Outside inpolygon','Outside inpoly2',Location='northwest')
title('P3')
subplot(2,2,4)
plot(xy(:, 1),xy(:, 2),'k.')
hold on
plot(p4,'FaceColor','none','EdgeColor','c', 'LineWidth',2)
plot(xy(Idx_Out_4, 1),xy(Idx_Out_4, 2),'xr')
title('P4')
plot(xy(Idx_Out_44, 1),xy(Idx_Out_44, 2),'ob')
legend('data','p4 Poly','Outside inpolygon','Outside inpoly2',Location='northwest')
% Method 1 Run Time = 5.3019
% Method 2 Run Time = 0.011318
As depicted in the figure above, the first subplot (P1) has some blue circles (inpoly2) falling within the boundaries, although these are supposed to be 'outside' points. Even worse is P3, where a huge number of datapoints flagged as 'outside' are indeed within the boundaries. Please note, taking P3 as an example, I am not interested in points falling outside/inside the smaller inner polygon. I only need those falling outside the major polygon.
Would you be able to help me finding a solution to correctly identify the points falling outside the boundaries?
Is there an alternative fast approach to find the 'outsiders'?
Any help would be greatly appreciated.
Thanks in advance!

请先登录，再回答此问题。

采纳的回答

Bruno Luong
Bruno Luong 2023-10-2
编辑：Bruno Luong 2023-10-3
The reason that inpoly2/insidepoly fail because the polygons p1 and p3 have holes, therefore px.Vertices have NaN values.
Here the fix: https://www.mathworks.com/matlabcentral/fileexchange/27840-2d-polygon-interior-detection?s_tid=srchtitle
EDIT1: simplified inpolwithhole
EDIT2: fix a bug
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Method 1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tic
is_outside1 = ~inpolygon(xy(:, 1), xy(:, 2), p1.Vertices(:, 1), p1.Vertices(:, 2));
Idx_Out_1 = find(is_outside1);
is_outside2 = ~inpolygon(xy(:, 1), xy(:, 2), p2.Vertices(:, 1), p2.Vertices(:, 2));
Idx_Out_2 = find(is_outside2);
is_outside3 = ~inpolygon(xy(:, 1), xy(:, 2), p3.Vertices(:, 1), p3.Vertices(:, 2));
Idx_Out_3 = find(is_outside3);
is_outside4 = ~inpolygon(xy(:, 1), xy(:, 2), p4.Vertices(:, 1), p4.Vertices(:, 2));
Idx_Out_4 = find(is_outside4);
toc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Method 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tic
[is_outside11] = inpolwithhole(xy, p1.Vertices);
Idx_Out_11 = find(is_outside11==0);
[is_outside22] = inpolwithhole(xy, p2.Vertices);
Idx_Out_22 = find(is_outside22==0);
[is_outside33] = inpolwithhole(xy, p3.Vertices);
Idx_Out_33 = find(is_outside33==0);
[is_outside44] = inpolwithhole(xy, p4.Vertices);
Idx_Out_44 = find(is_outside44==0);
toc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Plot %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure
subplot(2,2,1)
plot(xy(:, 1),xy(:, 2),'.','Color',0.8+[0 0 0])
hold on
plot(p1,'FaceColor','none','EdgeColor','r', 'LineWidth',2)
plot(xy(Idx_Out_1, 1),xy(Idx_Out_1, 2),'xr')
plot(xy(Idx_Out_11, 1),xy(Idx_Out_11, 2),'ob')
legend('data','p1 Poly','Outside inpolygon','Outside insidepoly',Location='northwest')
title('P1')
subplot(2,2,2)
plot(xy(:, 1),xy(:, 2),'.','Color',0.8+[0 0 0])
hold on
plot(p2,'FaceColor','none','EdgeColor','m', 'LineWidth',2)
plot(xy(Idx_Out_2, 1),xy(Idx_Out_2, 2),'xr')
plot(xy(Idx_Out_22, 1),xy(Idx_Out_22, 2),'ob')
legend('data','p2 Poly','Outside inpolygon','Outside insidepoly',Location='northwest')
title('P2')
subplot(2,2,3)
plot(xy(:, 1),xy(:, 2),'.','Color',0.8+[0 0 0])
hold on
plot(p3,'FaceColor','none','EdgeColor','g', 'LineWidth',2)
plot(xy(Idx_Out_3, 1),xy(Idx_Out_3, 2),'xr')
plot(xy(Idx_Out_33, 1),xy(Idx_Out_33, 2),'ob')
legend('data','p3 Poly','Outside inpolygon','Outside insidepoly',Location='northwest')
title('P3')
subplot(2,2,4)
plot(xy(:, 1),xy(:, 2),'.','Color',0.8+[0 0 0])
hold on
plot(p4,'FaceColor','none','EdgeColor','c', 'LineWidth',2)
plot(xy(Idx_Out_4, 1),xy(Idx_Out_4, 2),'xr')
title('P4')
plot(xy(Idx_Out_44, 1),xy(Idx_Out_44, 2),'ob')
legend('data','p4 Poly','Outside inpolygon','Outside insidepoly',Location='northwest')
%%
function b = inpolwithhole(xy, P)
isrownan = any(isnan(P),2);
if any(isrownan)
np = size(P,1);
rows = (1:np+1)';
rownan = find(isrownan);
rownannext = rownan+1;
% wrap each component back to the first element
rows([rownan; end]) = [1; rownannext];
% insert 1s afer each component
d = cumsum(accumarray(rownannext,2,[np+1,1],[],1));
rexpand = accumarray(d,rows,[],[],1);
Pe = P(rexpand,:);
else
Pe = P;
end
b = insidepoly(xy, Pe);
end
Timing on my PC for inpolygon and insidepoly (mex compiled version) are respectively:
  • Elapsed time is 1.689946 seconds.
  • Elapsed time is 0.004488 seconds.
  3 个评论
显示 1更早的评论
Bruno Luong
Bruno Luong 2023-10-5
Elapsed time is 0.006547 seconds.
for this code
function b = inpoly2holeremove(xy, p)
p=rmholes(p);
[is_outside11, ~] = inpoly2(xy, p.Vertices);
b = find(is_outside11==0);
end
Bruno Luong
Bruno Luong 2023-10-5
If you don't care about holes, this also work, still fast
Elapsed time is 0.004306 seconds.
function b = insidepolyholeremove(xy, p)
b = insidepoly(xy, rmholes(p).Vertices);
end

请先登录，再进行评论。

更多回答（1 个）

Steven Lord
Steven Lord 2023-10-5
How are you representing your polygon? If you've created them as polyshape objects try creating the polygon once as a polyshape object then calling isinterior for each of your data sets on the polyshape.
  1 个评论
Bruno Luong
Bruno Luong 2023-10-5
移动：Bruno Luong 2023-10-5
Yuck, no. Performance wise it is absurdly slow. I just test and isinterior is even slower than inpolygon used by Simon; it takes
Elapsed time is 3.789734 seconds.
600-900 tile slower than our non-official code.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Computational Geometry 的更多信息

标签

产品


版本

R2022b

Community Treasure Hunt

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

Start Hunting!

Translated by