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

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]

You might try mine https://www.mathworks.com/matlabcentral/fileexchange/27840-2d-polygon-interior-detection?s_tid=srchtitle

[Simone A.]

Hi @Bruno Luong Thanks for gettin back and thanks for providing another grat alternative tool!

I have already tried it but inpoly2, overall, seems more consistent with inpolygon.

Answer 1

[Bruno Luong]

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.

[Simone A.]

Hi @Bruno Luong Thanks a lot for your help! That indeed fixed the issue, great! Another approach I found (which might be slower however; would be great if you could test it on you (faster) machine, if you have some spare time, just for reference) is the one below :

p1=rmholes(p1);
[is_outside11, ~] = inpoly2(xy, p1.Vertices);
Idx_Out_11 = find(is_outside11==0);
% [Rest of the code...]

Once again, thanks a lot for your time and patience!

[Bruno Luong]

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]

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

Answer 2

[Steven Lord]

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.

[Bruno Luong]

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.