点和多边形关系判定算法

点和多边形关系判定是一个很常用的问题。常用算法是射线法和夹角法。

射线法=从该点假想一条通向无穷远的射线，通过判断射线与多边形交点的个数的奇偶来确定点与多边形的位置关系。

夹角法=计算多边形任意两相邻顶点对待求点的夹角，通过判断夹角和为0, 180, 360来确定点与多边形的关系。

由于计算夹角需要使用反三角函数，效率较低，所以国内的论文通常的都是倾向于射线法的。射线法的问题在于，所做的射线可能与多边形某一边重叠。判断一条射线与给定线段是否有交，似乎也不是那么简单。

让我们看看matlab是怎么做的1

function [relation] = inpolygen(x, y, xv, yv)
% 判断点(x, y)与由xv, yv确定的多边形的关系

Nv = length(xv); % 多边形顶点数
xv = xv - x; % 根据x,y 平移坐标系
yv = yv - y; % 此后只需判断(0,0)是否在多边形中即可

% 判断多边形顶点所在的象限0, 1, 2, 3分别表示1-4象限
posX = xv > 0;
posY = yv > 0;
negX = ~posX;
negY = ~posY;
quad = (negX & posY) + 2*(negX & negY) + 3*(posX & negY); % 很巧妙

% 两组循环变量
m = 1:Nv-1;
mp1 = 2:Nv;

% 计算相邻顶点向量积的符号
signCrossProduct = sign(xv(m) .* yv(mp1) - xv(mp1) .* yv(m));
% 计算相邻顶点数量积
dotProduct = xv(m) .* xv(mp1) + yv(m) .* yv(mp1);

% 相邻顶点所在象限的差
diffQuad = diff(quad);
% 修复象限差，使其属于{-2,-1,0,1,2}
idx = (abs(diffQuad) == 3);
diffQuad(idx) = -diffQuad(idx)/3;
idx = (abs(diffQuad) == 2);
diffQuad(idx) = 2*signCrossProduct(idx);

% 叉积等于0表示两向量平行, 说明(0,0)经过这条边所在的直线
% 此时点积<=0证明(0,0)在此线段上
if any((signCrossProduct == 0) & (dotProduct <= 0));
    relation = 'on';
    return;
end

% 象限差的和不等于0，在多边形内
if (sum(diffQuad) ~= 0)
    relation = 'in';
    return;
end

relation = 'out';
return;

整个算法的实现非常棒，把它转换为c, java代码都很容易。使用象限差来代替夹角，也避免了计算反三角函数，效率应该很不错。

回头看看维普上的论文，唉，不由得让人感慨。教授和商业公司都是为财，怎么差距就这么大呢？我还是建议经常会接触到这种基础计算的朋友，手上保留一份matlab的m文件，在%matlab%\toolbox\matlab下。不在校园的，我推荐使用开源的R。阅读它们的算法，常常可以带给你惊喜。

Update Dec.19 2008 : 移植了一个cpp的实现，在原有matlab的基础上进一步做了一些优化。10,000,000边的多边形在我的电脑上只需要不到0.2秒。从Google过来的点这里下载。

注1：matlab的代码更长，能够计算许多点和许多多边形两两之间的关系，这里我将其略去，并做了一些适合阅读的修改。

Post by : cooper

Post under : Coding, , , , , ,

on December 17th, 2008 | 5 Comments »