如何简单计算图像和多边形之间的重叠区域？

Question

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

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

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

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

然而，只有当垃圾桶的中心点包含在多边形内时，它才会返回true，也就是说，它返回了下面图片中的红色形状。然而，我需要更接近绿色形状（尽管还不完整）的解决方案。

为了获得绿色图形，我执行了四次对inpolygon的调用。对于每个比较，我将点网格向NE、NW、SE或SW移动了1/2。这相当于测试垃圾桶的角落是否在多边形内。

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);

虽然这给我提供了部分解决方案，但在一个顶点包含在一个箱子中但没有一个箱子角落的情况下它会失败。

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

使用以下内容绘制了此图：

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

- slayton

---

不得不离开电脑，但想提供一个通用解决方案，可能会有所帮助。首先将网格放大到22的倍数，以定义密度等于或大于顶点使用的密度的区域-这将消除角落问题。然后要缩小到22乘22的网格，只需除以您放大的相同因子，将顶部/左侧的点向下取整，将底部/右侧的点向上取整即可。希望能有所帮助。 - Salain

Answer 1

建议使用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;

- klurie

Answer 2

轻微改进

首先，为了简化您的“部分解决方案”——您只是看了图形的角落。如果您考虑23x23个角落的网格而不是22x22个点的网格（它们将被偏移-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) 和 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（对于边缘y = 1.5）。因此，我们可以将(1, 1)和(1, 2)的像素标记为我们集合的一部分。
• x = 1.75 被四舍五入为2（对于边缘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)的像素。

完成！

好吧，有点。这将给您边缘，但它不会填充多边形的主体。但是，您可以将其与先前（红色）结果结合使用，以获得完整集合。

- cloudfeet

Answer 3

这个伪代码算法怎么样？

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中检查相交是否容易。

- Superbest

Answer 4

我建议在图像处理工具箱中使用poly2mask，它可以做你想要的事情，而且也基本上符合你和Salain所建议的。

- Juhl

---

是的，我已经尝试过poly2mask了。它给出的结果与我在问题中描述的inpolygon调用相同。 - slayton

Answer 5

首先，我为这个例子定义了一个低分辨率的圆形

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函数将我们的高分辨率网格中的点四舍五入到最近的低分辨率等价点。然后，我们通过使用'rows'限定符调用unique函数以向量样式或成对方式删除所有重复项。就这样。

要查看结果，请执行以下操作：

[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

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

- Salain

Answer 6

好吧，我想我来晚了，严格来说赏金时间是到明天的; )。但是这里是我的尝试。首先，一个标记包含/接触点的单元格的函数。给定间距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函数获得多边形的内部。如果您需要有关代码的更多详细信息，请告诉我。

- angainor