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

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

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

1. 构造一个以两个多边形顶点为元素的集合 S
2. 对于多边形 A 的每条边 e1
  1. 对于多边形 B 的每条边 e2
    1. 如果 e1 与 e2 存在交点
      1. 将交点放入集合 S 中
3. 判断集合 S 中的每个元素,若不在 A 或 B 之内,则从 S 中剔除

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

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

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

inpolygon 函数用法如下

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

其中 xv, yv 表示构成多边形各顶点的一组列向量,与 polyxpoly 函数中的 x, y 类似;X, Y 为需要判断的点的坐标构成的列向量,例如 \(X = {[1,3,5,7]^T}\),\(Y = {[8,6,4,2]^T}\),则代表需要判断是否位于多边形内的四点分别是\((1,8)\),\((3,6)\),\((5,4)\) 和 \((7,2)\)。返回的 IN 为一逻辑列向量,长度与 X、Y 相等,若为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坐标,均为列向量
% ********************************************************** %
% 1.Let S be the set of vertices from both polygons.
% 2.For each edge e1 in polygon 1
% -- 1.For each edge e2 in polygon 2
% ---- 1.If e1 intersects with e2
% ------ 1.Add the intersection point to S
% 3.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); % 上一步得到的顶点坐标有可能不是按逆时针顺序的,因此利用 convhull 函数对它们对凸包
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 坐标构成的一对列向量,即 \(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

效果如下。


1 Comment

  1. Bone 2014-12-25 Reply

    bbnhs,nzbzdwhcbn!!!哈哈哈!

Leave a reply

Your email address will not be published.