求多边形(Convex Hull)交集的 MATLAB 实现

这几天在看色域映射(Gammut Mapping)内容时碰到了求凸包交集(Intersection of Convex Hulls)的问题,有一个比较巧妙的算法可以轻松地解决这个问题,这里记录备用。实际上这个算法并不限制于凸多边形,对于存在内陷的多边形也能给出交集。将一些二维条件置换为三维参数,可以很容易地把这个算法扩展到三维空间中,即求多面体的交集(Intersection of Polyhedron)。

如上图所示,深蓝色和红色代表两个多边形 A 和 B,中间浅蓝色部分即是两个多边形的交集。在 MATLAB 中,我们已知多边形 A 和 B 各顶点的坐标,现在要求的是浅蓝色交集构成的多边形的顶点坐标。算法其实很简单,遵循以下流程:

构造一个以两个多边形顶点为元素的集合 S

for 多边形 A 的每条边 as e1
    for 多边形 B 的每条边 as e2
        if (e1 与 e2 存在交点)
            将交点放入集合 S 中

for 集合 S 中的每个元素 as e
    if (e 属于 A 且 e 属于 B)
        % 不做处理,即保留 e
    else
        将 e 从 S 中删除

返回 S

如果需要求并集(Union),只需要把第9行的 e 属于 A e 属于 B 改为 e 属于 A e 属于 B 即可。

这里涉及两个问题:如何求两条线段交点坐标,以及如何判断一点是否位于多边形之内。所幸的是,这两个问题 MATLAB 都提供了现成的函数。polyxpoly 函数(需安装绘图工具箱)能够计算多边形的交点坐标,这里我们仅仅用它来计算两条线段的交点坐标;而 inpolygon 函数能够判断一点是否位于多边形内。polyxpoly 函数用法如下

[xi, yi] = polyxpoly(x, y, X, Y);

其中 x 为第一个多边形各顶点 x 坐标构成的列向量,y 为第一个多边形各顶点 y 坐标构成的列向量,XY 同理。对于两条线段来说,$x = {[{x_1},{x_2}]^T}$,$y = {[{y_1},{y_2}]^T}$,$X = {[{X_1},{X_2}]^T}$,$Y = {[{Y_1},{Y_2}]^T}$,如下图所示。对于多边形来说 $x_i$、$y_i$ 是代表交点坐标的一对列向量,但是对于两条线段来说 $x_i$、$y_i$ 则仅仅表示交点的坐标。

inpolygon 函数用法如下

IN = inpolygon(X, Y, xv, yv);

其中 xvyv 表示构成多边形各顶点的一组列向量,与 polyxpoly 函数中的 xy 类似;XY 为需要判断的点的坐标构成的列向量,例如 $X = {[1,3,5,7]^T}$,$Y = {[8,6,4,2]^T}$,则代表需要判断是否位于多边形内的四点分别是$(1,8)$,$(3,6)$,$(5,4)$ 和 $(7,2)$。返回的 IN 为一逻辑列向量,长度与 XY 相等,若为1则表示对于的 $(X, Y)$ 位于多边形内(或恰好在多边形上),为0则表示在多边形外。例如先创建一个五边形:

t = linspace(0,2*pi,6);
xv = cos(t);
yv = sin(t);
xv = [xv,xv(1)];
yv = [yv,yv(1)];

然后随机产生100个坐标点

X = randn(100,1);
Y = randn(100,1);

判断这100个点是否位于五边形内,返回一列向量 IN,为1则代表是,为0代表否

IN = inpolygon(X, Y, xv, yv);

让100个点中位于五边形内的用红色星号表示,五边形外的用蓝色圆圈表示,绘制图像

plot(xv, yv, 'k', X(IN), Y(IN), '*r', X(~IN), Y(~IN), 'ob');

结果如下图所示

有了以上两个函数,就可以进行多边形求交集了。我一般习惯用 convhull 函数生成多边形,即随机生成若干个点,然后求这些点的最小凸包就得到了一个凸多边形。Qhull 算法给出了求凸包(二维至九维)的快速实现。假设已经得到了多边形 A 的顶点 x, y 坐标的一对列向量 poly1_x、poly1_y,以及多边形 B 的一对列向量 poly2_x、poly2_y,则它们的交集的各顶点坐标构成的一对列向量为 [ints_x,ints_y]。以下是完整代码。

function [ints_x, ints_y] = polygon_intersect(poly1_x, poly1_y, poly2_x, poly2_y)
% 求两个个凸多边形的交集
% poly1_x,poly1_y 分别为第一个多边形的各个顶点的x,y坐标,均为列向量
% poly2_x,poly2_y 分别为第二个多边形的各个顶点的x,y坐标,均为列向量
% ********************************************************** %
% Let S be the set of vertices from both polygons.
% For each edge e1 in polygon 1
%   For each edge e2 in polygon 2
%     If e1 intersects with e2
%       1. Add the intersection point to S
% Remove all vertices in S that are outside polygon 1 or 2
% ********************************************************** %
S(:,1) = [poly1_x; poly2_x]; % 将两个多边形的坐标存入 S 中,顺序无所谓
S(:,2) = [poly1_y; poly2_y];
num = size(poly1_x, 1) + size(poly2_x, 1) + 1;
for i = 1:size(poly1_x, 1) - 1
    for j =1:size(poly2_x, 1) - 1
        X1 = [poly1_x(i); poly1_x(i+1)];
        Y1 = [poly1_y(i); poly1_y(i+1)];
        X2 = [poly2_x(j); poly2_x(j+1)];
        Y2 = [poly2_y(j); poly2_y(j+1)];
        [intspoint_x, intspoint_y] = polyxpoly(X1, Y1, X2, Y2); % 求两条线段交点的x,y坐标
        if ~isempty(intspoint_x) % 若两条线段无交点则跳至下一组线段,若有交点则将交点的x,y坐标存至S中
            S(num, 1) = intspoint_x;
            S(num, 2) = intspoint_y;
            num = num + 1; % 存入 S 后往下递推一行
        end
    end
end
IN = inpolygon(S(:,1), S(:,2), poly1_x, poly1_y);
S(IN == 0, :) = []; % 剔除掉不位于多边形 A 中的顶点坐标
IN = inpolygon(S(:,1), S(:,2), poly2_x, poly2_y);
S(IN == 0, :) = []; % 剔除掉不位于多边形 B 中的顶点坐标
X = S(:, 1); % 得到交集多边形的各个顶点坐标
Y = S(:, 2);
k = convhull(X, Y);
ints_x = X(k);
ints_y = Y(k);
plot(poly1_x, poly1_y, 'r', poly2_x, poly2_y, 'b', ints_x, ints_y, 'k');

最终效果如下图所示,中间紫色的多边形即表示多边形 A 和 B 的交集。

如果需要求多个多边形的交集,则可以循环调用 polygon_intersect 函数。当某两个多边形之间交集为零时,可以考虑以同样倍数增大这两个多边形使得它们之间存在交集(下面的代码并没有考虑到这点)。这里我将输入设定为元胞数组,其第 i 项为第 i 个多边形的 x, y 坐标构成的一对列向量,即 $\mathrm{polygon}\{i\} = \left[ {{{[{x_1},{x_2},{x_3},…]}^T},{{[{y_1},{y_2},{y_3},…]}^T}} \right]$。

function [ints_x, ints_y] = Multipolyints(polygon)
if ~strcmp(class(polygon), 'cell') % 检查输入是否为元胞数组
    error('The input must be a cell.');
end
for i = 2:length(polygon) % 循环调用 polygon_intersect 函数
    polygon1 = polygon{1};
    polygon2 = polygon{i};
    [ints_x, ints_y] = polygon_intersect(polygon1(:,1), polygon1(:,2), polygon2(:,1), polygon2(:,2));
    polygon1 = [ints_x, ints_y];
end

效果如下。

About the author

Jueqin

本作品以 CC BY-NC-ND 许可协议进行发布。

如果您认为文章对您有用的话,不妨请我喝一杯咖啡?

12 comments

  • 用matlab自带的函数polyshape和intersect函数,参考https://ww2.mathworks.cn/help/matlab/ref/polyshape.intersect.html?s_tid=gn_loc_drop描述。精确度不错。

    • 是的,新版本下自带函数已经很好用了,不过写这篇文章的时候几何学工具箱还不怎么齐全。

  • S(num,1) = intspoint_x;
    如果两条线段完全重合,polyxpoly会返回两个交点坐标。
    此时,num报错

  • 有个问题,如果相交的两条线中有一条是平行于坐标轴的,比如y=3,交点坐标可能会出现y=3.000000000000001这种,大概会有10e-15数量级的余量,导致后面inpolygon函数判断的时候把这个交点删除掉,最后convhull的时候可能因为点不足3个而报错。不知道是我个人的问题,还是polyxpoly这个函数的问题。