直线是穿过顶点，还是切过顶点？

最新推荐文章于 2024-10-05 22:19:52 发布

多彩电脑

最新推荐文章于 2024-10-05 22:19:52 发布

阅读量1.2k

点赞数 20

分类专栏：有意思的数学知识文章标签：算法几何学 matlab

本文链接：https://blog.csdn.net/weixin_41640959/article/details/141749299

版权

有意思的数学知识专栏收录该内容

4 篇文章 0 订阅

订阅专栏

关于判断点是否在多边形内，有一篇博文[1]，详谈判断点在多边形内的七种方法（最全面） hdu1756 hrbust1429 为例_判断一个点在不在多边形内-CSDN博客

该博文介绍了多种判断点是否在多边形内的方法。其中该博文介绍的第一种方法，是作射线统计与多边形的交点数。该方法有特殊情况，例如该射线与多边形的某个顶点相交，是否应当计数。

注：该图由博文[1]给出

本文对于射线刚好和顶点相交的时该顶点是否应当被计数进行进一步分析，同时给出两种方法进行判断。与此同时，会附上判断算法的Matlab代码。

一、射线经过的顶点是否应被计数的判断几何原理

（一）一般情形

博文[1]的第一种方法，每一次的计数，实质上是该射线在向右射的过程中，内外状态的改变

以该图为例：

（1）对于以P点为起点的射线，假定某可动点从P出发，在蓝色a段上，该可动点仍然在多边形内；但当到达Q点时，情况出现了变化，可动点"碰"到了多边形，从在多边形内开始向在多边形外转化。可动点继续向右滑动，在红色b段上（从点Q开始向右无限延伸），可动点永远处于多边形外。显然，P在多边形内，转化了一次后成为多边形外。

（2）对于以G点为起点的射线，假定某可动点从G出发，在蓝色c段上，该可动点仍然在多边形外；但到了黄色d段上，该可动点处于多边形内；继续滑动，进入红色e段后，可动点永远处于多边形外。显然，G在多边形外，转化了一次后成为多边形内，再转化了一次后成为多边形外。

由于多边形是有限大的图形，任何射线上的可动点从起点出发，从某一点开始都会永远处于多边形外，因此最终状态必为"外"。所以，内外转化次数为奇数，则该点处于多边形内；内外转化次数为偶数（包括0），则该点处于多边形外。

（二）特殊情形

若射线不和顶点相交，则和边碰撞时会发生一次内外转化。但是，若直接和顶点相交，则内外转化未必发生。

如图，可动点在以M为起点的射线，在N点和多边形的边相交，内外转化发生了，从外到内。但在经过O顶点时，经过前在内部，经过后仍在内部，没有发生内外转化。显然，这一次不应计数。但是，可动点在以W为起点的射线，在经过U顶点时，发生了内外转化，从内到外。

由图可见，射线 $\vec{MO}$ 和射线 $\vec{WU}$ 的区别在于，线段OS和OE在射线 $\vec{MO}$ 的同侧，所以相当于以M为起点的射线在点O和折线SOE相切；而线段UD和UV在 $\vec{WU}$ 的异侧，所以相当于以W为起点的射线在点U穿过折线RUV。

因此，若以A为起点的射线和多边形的顶点B相交，假设和顶点B相连的顶点C和顶点D在射线A的同侧，即射线A和折线CBD在D相切，则不发生内外转化，故而在博文[1]介绍的第一个方法中，不应计数。但若C和D在射线A的异侧，即射线A在D穿过折线CBD，则发生内外转化，应计数。

对于更特殊的情况，若C或D在射线A上，则忽略该边，将点B与C或D（射线上）的邻点连接起来，再进行判断。

可用此方法将其转化为和顶点相切或穿过顶点的情况。

二、判断经过顶点的直线与折线是否相切或穿过

在第一章已经说明了在射线法[1]中判断射线经过顶点是否应当计数是根据射线是否和顶点及其两边的线段构成的折线相切还是穿过折线判断的。本章介绍两种判断经过折线段顶点的直线和该折线段的关系的方法。

如图，第一张图中，直线CD和折线ACB相切。第二种图中，直线CD穿过折线ACB。判断的方式是，如果点A和点B在直线CD的同侧，则直线CD和折线ACB相切；如果点A和点B在直线CD的异侧，则直线CD穿过折线ACB。

关于同侧还是异侧的判断，本文给出两种方式：

（一）解析几何法

在平面直角坐标系中，一条直线可以表达成ax+by+c=0的形式。而如果点 $A(x_A,y_A)$ 不在直线上，而在直线的一侧，且 $ax_A+by_A+c>0$ ，则如果点 $B(x_B,y_B)$ 在直线的另一侧，那么 $ax_B+by_B+c<0$ 。反之亦然。对于直线上的点 $C(x_C,y_C)$ ， $ax_C+by_c+c=0$ 。

因此，要判断点A和点B在直线CD的同侧或异侧，或在直线CD上，步骤如下：

（1）计算直线CD的解析表达式ax+by+c=0，其中 $a=y_D-y_C$ ， $b=-(x_D-x_C)$ ， $c=y_Cx_D-x_Cy_D$ [2]。

（2）计算 $\Delta=(ax_A+by_A+c)(ax_B+by_B+c)$ 。若 $\Delta>0$ ，则点A和点B在直线CD的同侧；若 $\Delta<0$ ，则点A和点B在直线CD的异侧； $\Delta=0$ ，则点A或点B在直线CD上。

实现该判断的功能Matlab代码如下：

（1）该代码是一个判断直线CD和折线ACB相切，还是穿过折线ACB的函数。若返回1则穿过，若返回0则相切，若返回-1表明A或B至少有一个在直线CD上。该函数在一开始会判断A，B，C，D是否为4个不同的点。

function isCrossing = crossingVia(A, B, C, D)
%A is [x_A, y_A]. Others also
%isCrossing: 1, crossing. 0, tangent. -1, overlap

%Precondition: ABCD do not overlap
points = [A;B;C;D];
for ii = 1:4
    for jj = (ii+1):4
        mat1 = points(ii,:);
        mat2 = points(jj,:);
        if isequal(mat1, mat2)
            disp(['node ', num2str(ii), ' and node ', num2str(jj), ' are same!']);
            return
        end
    end
end

%Here we get the line CD's equation line_a * x + line_b * y + line_c = 0
line_a = D(2) - C(2);
line_b = -(D(1) - C(1));
line_c = C(2) * D(1) - C(1) * D(2);

A_with_CD = line_a * A(1) + line_b * A(2) + line_c;
B_with_CD = line_a * B(1) + line_b * B(2) + line_c;

if A_with_CD * B_with_CD > 0 %same side
    isCrossing = 0;
elseif A_with_CD * B_with_CD < 0 %different side
    isCrossing = 1;
else %overlap
    isCrossing = -1;
end
end

（2）接下来，为了直观地展示线段、直线的位置关系，用Matlab进行作图。

线段作图比较简单，用Matlab里的plot函数

function drawLine(pointA, pointB)
%This function draws a segment between pointA and pointB
plot([pointA(1) pointB(1)], [pointA(2) pointB(2)]);
end

（3）直线是向两边无限延长的线段，故需要延长线段，根据图的大小延长。博客[3]提供了示例代码及其说明。这里按照该博客提供的思路，定义一个函数画延长的直线。函数里，xlim是直线在x方向的范围，ylim是y方向的范围。注意程序里的第二个if语句，因为Matlab函数polyfit不适用于竖直线（斜率为 $\infty$ ），所以对于竖直线，需另外写程序。

function drawExtendedLine(pointA, pointB, xlim, ylim)
%This function draws a line passing pointA and pointB
%xlim, ylim is the x, y range of the line. From small to big
if isequal(pointA, pointB)
    disp("Node A is same as node B")
    return;
elseif (xlim(1) >= xlim(2)) || (ylim(1) >= ylim(2))
    disp("range must be from small to big")
    return;
end
if pointA(1) ~= pointB(1)
    polyRes = polyfit([pointA(1), pointB(1)], [pointA(2), pointB(2)], 1);
    k = polyRes(1);%(pointB(2) - pointA(2)) / (pointB(1) - pointA(1));
    b = polyRes(2);%pointA(2) - k * pointA(1);
    leftPoint = [xlim(1), k * xlim(1) + b];
    rightPoint = [xlim(2), k * xlim(2) + b]; %get the equation of the line
else
    leftPoint = [pointA(1) ylim(1)];
    rightPoint = [pointA(1) ylim(2)]; %vertical line (slope = inf)
end
drawLine(leftPoint, rightPoint)
end

（4）另外，在图上可以标注点的名称，用Matlab函数text。

function labelNode(point, name)
text(point(1),point(2),name)
end

（5）最后，定义A，B，C，D四个点，并用Matlab把点和线都画出来，另外用解析几何判断A，B在直线CD的同侧还是异侧。

hold off %clear graph
nodeA = [0,5];
nodeB = [1,9];
nodeC = [5,7];
nodeD = [8,5];%position of four nodes
drawLine(nodeA,nodeC)
xlim([0,10])
ylim([0,10])
hold on
drawLine(nodeC,nodeB)
drawExtendedLine(nodeC,nodeD,[0,10],[0,10])
%Now label each node
labelNode(nodeA, 'A')
labelNode(nodeB, 'B')
labelNode(nodeC, 'C')
labelNode(nodeD, 'D')
%Finally, check whether the line CD crosses poly line ACB
crossingOrNot = crossingVia(nodeA, nodeB, nodeC, nodeD);
if crossingOrNot == 1
    disp("CD crosses polyline ACB")
elseif crossingOrNot == 0
    disp("CD is tangent to polyline ACB")
elseif crossingOrNot == -1
    disp("CD is overlapping with either AC or BC")
else
    disp("Unknown error")
end

结果如下：

CD和折线ACB相切。

（二）向量法

两个起点相同的向量的叉积（即向量积）可用于判断两个向量之间的位置关系[4]。叉积的定义，以及运算方式，见[4]。简单地说，对于平面向量 $(x_1,y_1)$ 和 $(x_2,y_2)$ ， $(x_1,y_1)\times (x_2,y_2)=(x_1y_2-x_2y_1)\vec{k}$ ，其中 $\vec{k}$ 是z轴的单位向量。令 $\Delta=x_1y_2-x_2y_1$ ，则当 $\Delta>0$ 时，向量 $(x_2,y_2)$ 在 $(x_1,y_1)$ 的逆时针方向；当 $\Delta<0$ 时，向量 $(x_2,y_2)$ 在 $(x_1,y_1)$ 的顺时针方向；当 $\Delta=0$ 时， $(x_2,y_2)$ 与 $(x_1,y_1)$ 在同一直线上[4]。

以C为中心，作向量 $\vec{CA}$ ， $\vec{CB}$ ， $\vec{CD}$ ：

如图可知，当 $\vec{CD}$ 在 $\vec{CA}$ ， $\vec{CB}$ 的同一侧时，直线CD和折线ACB相切，此时 $\Delta_{\vec{CA}\times\vec{CD}}* \Delta_{\vec{CB}\times\vec{CD}} >0$ 。而当 $\vec{CD}$ 在 $\vec{CA}$ ， $\vec{CB}$ 的不同侧时，直线CD穿过折线ACB，此时 $\Delta_{\vec{CA}\times\vec{CD}}* \Delta_{\vec{CB}\times\vec{CD}}<0$ 。若重合，则 $\Delta_{\vec{CA}\times\vec{CD}}* \Delta_{\vec{CB}\times\vec{CD}}=0$ 。

$\vec{CD}=(x_D-x_C,y_D-y_C)$ ，其它向量类似。

实现该判断的功能Matlab代码如下：

（1）Matlab有一个函数cross可以计算叉积，计算时需输入2个三维向量。对于平面向量而言，第三维=0。计算的结果也是三位向量，对于平面向量而言，第三维的值即 $\Delta$ 。所以定义一个函数计算两个向量的叉积的 $\Delta$ 。

function scale = crossProductScale(vectorA, vectorB)
%This function calculate the scale part of the crossproduct of two 2D
%vectors (vectorA and vectorB) based on the assumption that they share the same origin. So if vectorB is on the CCW side of vectorA
%then scale>0, if on CW side then scale<0, and if in the same line then
%scale=0
crossVector = cross([vectorA 0], [vectorB 0]);
scale = crossVector(3);
end

（2）判断是否相切或穿过。首先根据点坐标，计算 $\vec{CA}$ ， $\vec{CB}$ ， $\vec{CD}$ 。然后计算叉积， $\Delta$ ，最后做出判断。

function isCrossingOrNot = checkWhetherCrossing(A, B, C, D)
%get the vectors first
vectorCA = A - C;
vectorCB = B - C;
vectorCD = D - C;
sameSide = crossProductScale(vectorCA, vectorCD) * crossProductScale(vectorCB, vectorCD);
if sameSide > 0
    isCrossingOrNot = 0;
elseif sameSide < 0
    isCrossingOrNot = 1;
else
    isCrossingOrNot = -1;
end
end

（3）画线段，直线，标注点的名称的代码前文已提到。

（4）Matlab有画箭头的程序，但要将箭头画在坐标系的指定位置，需自行定义函数进行转换。博文[5]对此提供了一种解决方案。在这里采用该方案定义画箭头函数。

function drawArrow(pointA, pointB)
%https://blog.csdn.net/u012366767/article/details/99568619
posAxes = get(gca, 'Position');
posX = posAxes(1);
posY = posAxes(2);
width = posAxes(3);
height = posAxes(4);
limX = get(gca, 'Xlim');
limY = get(gca, 'Ylim');
minX = limX(1);
maxX = limX(2);
minY = limY(1);
maxY = limY(2);
xForTwo = posX + [(pointA(1)-minX)/(maxX-minX) (pointB(1)-minX)/(maxX-minX)] * width;
yForTwo = posY + [(pointA(2)-minY)/(maxY-minY) (pointB(2)-minY)/(maxY-minY)] * height;
annotation("arrow", xForTwo, yForTwo)
end

（5）最后，定义A，B，C，D四个点，并用Matlab把点和线都画出来，另外用向量法判断A，B在直线CD的同侧还是异侧。

nodeA = [6,5];
nodeB = [7,5];
nodeC = [5,1];
nodeD = [4,7];
nodeCrossing = checkWhetherCrossing(nodeA, nodeB, nodeC, nodeD);
hold off
range = [0,10];
xlim(range)
ylim(range)
drawArrow(nodeC, nodeA)
hold on
drawArrow(nodeC, nodeB)
drawExtendedLine(nodeC, nodeD, range, range)
labelNode(nodeA, 'A')
labelNode(nodeB, 'B')
labelNode(nodeC, 'C')
labelNode(nodeD, 'D')
if nodeCrossing == 1
    disp("CD cross ACB")
elseif nodeCrossing == 0
    disp("CD tangent via ACB")
elseif nodeCrossing == -1
    disp("CD overlaps with either CA or CB")
else
    disp("Unknown error")
end

结果如下：