Help finding datapoints in x,y falling outside the boundaries of a polygon. Incongruencies between inpolygon and inpoly2, time issues, and possible solutions
6 views (last 30 days)
Show older comments
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!
2 Comments
Accepted Answer
Bruno Luong
on 2 Oct 2023
Edited: Bruno Luong
on 3 Oct 2023
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 Comments
Bruno Luong
on 5 Oct 2023
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
on 5 Oct 2023
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
More Answers (1)
Steven Lord
on 5 Oct 2023
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 Comment
Bruno Luong
on 5 Oct 2023
Moved: Bruno Luong
on 5 Oct 2023
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.
See Also
Categories
Find more on Computational Geometry in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!