7

我有一个封闭的非自相交多边形。它的顶点保存在两个向量 X 和 Y 中。最后 X 和 Y 的值被限制在 0 和 22 之间。

我想构造一个大小为 22x22 的矩阵，如果多边形的一部分与该 bin 重叠，则将每个 bin 的值设置为 true，否则为 false。

我最初的想法是生成一个由定义的点组成的网格[a, b] = meshgrid(1:22)，然后用来inpolygon确定网格的哪些点在多边形中。

[a b] = meshgrid(1:22);
inPoly1 = inpolygon(a,b,X,Y);

然而，只有当 bin 的中心包含在多边形中时，才会返回 true，即它返回下图中的红色形状。然而，更需要沿着绿色形状的线条（尽管它仍然是一个不完整的解决方案）。

为了获得绿色斑点，我对inpolygon. 对于每次比较，我将 NE、NW、SE 或 SW 点的网格移动了 1/2。这相当于测试 bin 的角是否在多边形中。

inPoly2 = inpolygon(a-.5,b-.5,X,Y) | inpolygon(a+.5,b-.5,X,Y) | inpolygon(a-.5,b+5,X,Y) | inpolygon(a+.5,b+.5,X,Y);

虽然这确实为我提供了部分解决方案，但在顶点包含在 bin 中但没有 bin 角的情况下它会失败。

有没有更直接的方法来解决这个问题，最好是产生更具可读性的代码的解决方案？

在此处输入图像描述

这个情节是用：

imagesc(inPoly1 + inPoly2); hold on;
line(a, b, 'w.');
line(X, Y, 'y);
4

6 回答 6

5

一种建议是使用 polybool 函数（在 2008b 或更早版本中不可用）。它找到两个多边形的交点并返回结果顶点（如果不存在顶点，则返回一个空向量）。为了在这里使用它，我们迭代（使用arrayfun）检查网格中的所有正方形，以查看polybool 的输出参数是否为空（例如，没有重叠）。

N=22;
sqX = repmat([1:N]',1,N);
sqX = sqX(:);
sqY = repmat(1:N,N,1);
sqY = sqY(:);

intersects = arrayfun((@(xs,ys) ...
      (~isempty(polybool('intersection',X,Y,[xs-1 xs-1 xs xs],[ys-1 ys ys ys-1])))),...
      sqX,sqY);

intersects = reshape(intersects,22,22);

这是生成的图像：

在此处输入图像描述

绘图代码：

imagesc(.5:1:N-.5,.5:1:N-.5,intersects');
hold on;
plot(X,Y,'w');
for x = 1:N
    plot([0 N],[x x],'-k');
    plot([x x],[0 N],'-k');
end
hold off;
于 2012-09-04T19:07:29.463 回答
1

略有改善

首先，为了简化你的“部分解决方案”——你所做的只是看着角落。如果不考虑 22x22 的点网格，您可以考虑 23x23 的角网格（它将与较小的网格偏移 (-0.5, -0.5)。一旦有了，您就可以在 22x22 网格上标记点在多边形中至少有一个角。

完整解决方案：

但是，您真正要查找的是多边形是否与每个像素周围的 1x1 框相交。这不一定包括任何角，但它确实要求多边形与盒子的四个边之一相交。

您可以使用以下算法找到多边形与包含框相交的像素的一种方法：

For each pair of adjacent points in the polygon, calling them pA and pB:
    Calculate rounded Y-values: Round(pA.y) and Round(pB.y)
    For each horizontal pixel edge between these two values:
        * Solve the simple linear equation to find out at what X-coordinate
          the line between pA and pB crosses this edge
        * Round the X-coordinate
        * Use the rounded X-coordinate to mark the pixels above and below
          where it crosses the edge
    Do a similar thing for the other axis

因此，例如，假设我们正在查看pA = (1, 1)and pB = (2, 3)。

  • 首先，我们计算了四舍五入的 Y 值：1 和 3。
  • 然后，我们查看这些值之间的像素边缘：y = 1.5和y = 2.5（像素边缘从像素偏移一半
  • 对于其中的每一个，我们求解线性方程以找到pA->pB与我们的边缘相交的位置。这给了我们：x = 1.25, y = 1.5和x = 1.75, y = 2.5。
  • 对于这些交点中的每一个，我们取舍入的 X 值，并用它来标记边缘任一侧的像素。
    • x = 1.25舍入为 1（对于 edge y = 1.5）。因此，我们可以将像素标记为我们集合(1, 1)的(1, 2)一部分。
    • x = 1.75舍入为 2（对于 edge y = 2.5）。因此，我们可以标记 和 处的(2, 2)像素(2, 3)。

这就是水平边缘的处理。接下来，让我们看看垂直的：

  • 首先我们计算四舍五入的 X 值：1 和 2
  • 然后，我们查看像素边缘。在这里，只有一个：x = 1.5。
  • 对于这条边，我们找到它与线pA->相交的地方pB。这给了我们x = 1.5, y = 2.
  • 对于这个交叉点，我们取舍入的 Y 值，并用它来标记边缘任一侧的像素：
    • y = 2舍入为 2。因此我们可以标记 和 处的(1, 2)像素(2, 2)。

完毕！

嗯，有点。这将为您提供边缘，但不会填充多边形的主体。但是，您可以将这些与您之前的（红色）结果结合起来以获得完整的集合。

于 2012-09-05T13:32:49.593 回答
1

这个伪代码算法怎么样：

For each pair of points p1=p(i), p2=p(i+1), i = 1..n-1
    Find the line passing through p1 and p2
    Find every tile this line intersects // See note
    Add intersecting tiles to the list of contained tiles

Find the red area using the centers of each tile, and add these to the list of contained tiles

注意：这条线需要一点点努力才能实现，但我认为有一个相当简单、众所周知的算法。

另外，如果我使用的是 .NET，我只需定义一个与每个网格图块对应的矩形，然后查看哪些与多边形相交。但是，我不知道在 Matlab 中检查交叉点是否容易。

于 2012-08-31T23:36:39.910 回答
1

我建议poly2mask在图像处理工具箱中使用，我认为它或多或少可以满足您的需求，并且或多或少可以满足您和 Salain 的建议。

于 2012-09-05T12:52:31.253 回答
0

好吧，我想我迟到了，虽然严格来说赏金时间是到明天；）。但这是我的尝试。首先，一个标记包含/触摸一个点的单元格的函数。给定一个间距为 lx, ly 的结构化网格和一组坐标为 (Xp, Yp) 的点，集合包含单元格：

function cells = mark_cells(lx, ly, Xp, Yp, cells)

% Find cell numbers to which points belong.
% Search by subtracting point coordinates from
% grid coordinates and observing the sign of the result.
% Points lying on edges/grid points are assumed
% to belong to all surrounding cells.

sx=sign(bsxfun(@minus, lx, Xp'));
sy=sign(bsxfun(@minus, ly, Yp'));
cx=diff(sx, 1, 2);
cy=diff(sy, 1, 2);

% for every point, mark the surrounding cells
for i=1:size(cy, 1)
    cells(find(cx(i,:)), find(cy(i,:)))=1;
end
end

现在，剩下的代码。对于多边形中的每个线段（您必须一个接一个地穿过这些线段），将该线段与网格线相交。使用给定的网格点坐标，分别针对水平线和垂直线小心地进行交叉，以避免数值不准确。对于找到的交点，我调用 mark_cells 将周围的单元格标记为 1：

% example grid
nx=21;
ny=51;
lx = linspace(0, 1, nx);
ly = linspace(0, 1, ny);
dx=1/(nx-1);
dy=1/(ny-1);
cells = zeros(nx-1, ny-1);

% for every line in the polygon...
% Xp and Yp contain start-end points of a single segment
Xp = [0.15 0.61];
Yp = [0.1 0.78];

% line equation
slope = diff(Yp)/diff(Xp);
inter = Yp(1) - (slope*Xp(1));

if isinf(slope)
    % SPECIAL CASE: vertical polygon segments
    % intersect horizontal grid lines
    ymax = 1+floor(max(Yp)/dy);
    ymin = 1+ceil(min(Yp)/dy);
    x=repmat(Xp(1), 1, ymax-ymin+1);
    y=ly(ymin:ymax);
    cells = mark_cells(lx, ly, x, y, cells);
else
    % SPECIAL CASE: not horizontal polygon segments
    if slope ~= 0
        % intersect horizontal grid lines
        ymax = 1+floor(max(Yp)/dy);
        ymin = 1+ceil(min(Yp)/dy);
        xmax = (ly(ymax)-inter)/slope;
        xmin = (ly(ymin)-inter)/slope;
        % interpolate in x...
        x=linspace(xmin, xmax, ymax-ymin+1);
        % use exact grid point y-coordinates!
        y=ly(ymin:ymax); 
        cells = mark_cells(lx, ly, x, y, cells);
    end

    % intersect vertical grid lines
    xmax = 1+floor(max(Xp)/dx);
    xmin = 1+ceil(min(Xp)/dx);
    % interpolate in y...
    ymax = inter+slope*lx(xmax);
    ymin = inter+slope*lx(xmin);
    % use exact grid point x-coordinates!
    x=lx(xmin:xmax); 
    y=linspace(ymin, ymax, xmax-xmin+1);
    cells = mark_cells(lx, ly, x, y, cells);
end

示例网格/段的输出： 输出

遍历所有多边形段会为您提供多边形“光环”。最后，使用标准的 inpolygon 函数获得多边形的内部。如果您需要有关代码的更多详细信息，请告诉我。

于 2012-09-10T22:12:42.673 回答
0

首先，我为此示例定义了一个低分辨率圆

X=11+cos(linspace(0,2*pi,10))*5;
Y=11+sin(linspace(0,2.01*pi,10))*5;

就像您的示例一样，它适合约 22 个单位的网格。然后，按照您的指示，我们声明一个网格网格并检查点是否在多边形中。

stepSize=0.1;
[a b] = meshgrid(1:stepSize:22);
inPoly1 = inpolygon(a,b,X,Y);

唯一不同的是，在您的原始解决方案采取一步的地方，这个网格可以采取更小的步骤。最后，在正方形的“边缘”内包含任何内容

inPolyFull=unique( round([a(inPoly1) b(inPoly1)]) ,'rows');

简单地采用我们的round高分辨率网格并将点适当地四舍五入到最接近的低分辨率等效值。unique然后，我们通过使用'rows'限定符调用以向量样式或成对方式删除所有重复项。就是这样

要查看结果，

[aOrig bOrig] = meshgrid(1:22);
imagesc(1:stepSize:22,1:stepSize:22,inPoly1); hold on;
plot(X,Y,'y');
plot(aOrig,bOrig,'k.');
plot(inPolyFull(:,1),inPolyFull(:,2),'w.'); hold off;

示例多边形

更改stepSize具有以速度和内存为代价改善结果的预期效果。

如果您需要结果与示例中的 inPoly2 格式相同，则可以使用

inPoly2=zeros(22);
inPoly2(inPolyFull(:,1),inPolyFull(:,2))=1

希望有帮助。我可以想到一些其他的方法来解决它，但这似乎是最简单的。

于 2012-09-04T15:12:51.837 回答