3

When determining points within polygons using MATLAB's inpolygon function, I find that the results are exactly correct for polygons drawn on linear axes but only approximately correct for polygons drawn on log-scale axes. Although my suspicions lean in favor of a MATLAB bug, it's possible I've overlooked something.

The following code reproduces the issue I have been experiencing with other data. The results are shown in the following image (the bottom set of panels are zoomed views of the top panels). One can appreciate that there are unlabeled points inside the polygon and labeled points outside the polygon, neither of which should occur, in the case of a polygon drawn on log-scale axes (right). In contrast, the polygon test is exact for polygons drawn on linear axes (left).

n=2E4;
x(:,1)=rand(n,1); y(:,1)=rand(n,1);
x(:,2)=lognrnd(0.5,0.25,n,1); y(:,2)=lognrnd(0.5,0.25,n,1);
for m=1:2
    subplot(1,2,m);
    scatter(x(:,m),y(:,m),'.'); hold on;
    if(m==2)
        set(gca,'xscale','log'); set(gca,'yscale','log');
    end
    p=impoly(gca);
    pc=getPosition(p);
    in=inpolygon(x(:,m),y(:,m),pc(:,1),pc(:,2));
    scatter(x(in,m),y(in,m),20);
end

enter image description here

Adriaan
  • 17,741
  • 7
  • 42
  • 75
user001
  • 1,850
  • 4
  • 27
  • 42

1 Answers1

6

I think you missed something: A line in normal scale is not a line in log scale. Your polygons are not properly drawn in the log scale, as you draw 2 points and put them together with a straight line.

Look at the real polygon in log space:

close all
clear

n=2e4;
x(:,1)=rand(n,1); y(:,1)=rand(n,1);
x(:,2)=lognrnd(0.5,0.25,n,1); y(:,2)=lognrnd(0.5,0.25,n,1);
for m=1:2
    subplot(1,2,m);
    scatter(x(:,m),y(:,m),'.'); hold on;
    if(m==2)
        set(gca,'xscale','log'); set(gca,'yscale','log');
    end
    p=impoly(gca);
    pc=getPosition(p);
    % plot polygon
    hold on
    for ii=1:size(pc,1)-1
        plot(linspace(pc(ii,1),pc(ii+1,1),100),linspace(pc(ii,2),pc(ii+1,2),100),'g')
    end
    plot(linspace(pc(end,1),pc(1,1),100),linspace(pc(end,2),pc(1,2),100),'g')
    in=inpolygon(x(:,m),y(:,m),pc(:,1),pc(:,2));
    scatter(x(in,m),y(in,m),20);
end

Look at this zoomed in result (click to enlarge):

enter image description here

This happens because the polygon is defined in euclidean space, and it is defined as points linked by lines. If you want to work in log space, things may get complicated. One way to numerically approximate it is the inverse of what I did for plotting. Create dense enough sampled straight line on log space, convert it to linear space, and define a high vertex polygon with the resulting points. Then use inpolygon.

Ander Biguri
  • 35,140
  • 11
  • 74
  • 120
  • Oh, that's a good observation! I overlooked it completely. When I run your code (on R2014b), the polygons don't superimpose. The polygon drawn in linear space is copied over to log-space but since the axes are different, it winds up being in the lower left corner of the log-space axes. – user001 Dec 19 '16 at 13:18
  • @user001 which polygons dont superimpose? I mean, the whole point of the code is showing that they dont in the second figure, because `log` space. I made the code and plot in 2014b – Ander Biguri Dec 19 '16 at 13:19