matlab 求向量的交集_求多边形（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 坐标构成的列向量，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}$，如下图所示。对于多边形来说 $x_i$、$y_i$ 是代表交点坐标的一对列向量，但是对于两条线段来说 $x_i$、$y_i$ 则仅仅表示交点的坐标。

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坐标，均为列向量

% ********************************************************** %

% 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

效果如下。