引言
例如说在用程序实现一款游戏时，有时候需要判断两个物体是不是发生了碰撞，如果两个物体没发生碰撞，那么一个物体的边上的任意点必定在另一个物体之外，这就可以抽象为一般的问题：给定一点$p_x$，给定一组确定一个$n$边形的端点集$p_1,p_2,\cdots,p_n$，判断点$p_x$是否位于多边形$p_1,\cdots,p_n$所围成的范围内.
点是否位于多边形内
如上图，图中的px是在多边形内部的，而px一撇显然是在多边形外部．像是「一个点是否在一个三角形内部」，「一个点是否在一个矩形内部」这样的问题，都是该问题的特殊情．为了简化问题，我们暂时限定所有点是位于同样一个平面，即$p_x, p_1, \cdots, p_n \in \mathbb{R}^{2}$，3维的情形可能稍微要复杂一些.
射线法
想象点是在一条直线轨道内向某个方向移动，点每穿过一条直线，计数器加1，那么，假如说点是在外部，这个点跨过边界时，要穿过一条线，这个点再走出边界时，还要再跨过一条线，因此，如果点是在外部，那么，从这个点做一条射线，这条射线与多边形所有边的交点的个数，必定是偶数个，因为，所有这些交点的产生原因可分为两种：点进入多边形穿过的边界，点走出多边形穿出的边界，点在外部，因此点进入的次数必定等于穿出的次数，因此所有交点的个数是偶数个；而如果，点是在多边形内部，那么，点如果要穿出多边形，它越过的边的个数肯定是奇数个.
射线法的实现
我们只以三角形的情形为例，演示如何通过程序来判断点是否在三角形内部，矩形和更一般的n边形，留给读者自己去尝试，首先呢，程序也是通过坐标来描述点的，如上图，有A,B,C,D,E这5个定点，另外一个动点不妨记做F，坐标分别记做$(x_A,y_A),(x_B,y_B),(x_C,y_C),(x_D,y_D),(x_E,y_E),(x_F,y_F)$.
首先，从出发点D做一条90度射线，这条射线上的所有点$(x,y)$必定满足
$$
\begin{bmatrix}
x \\
y
\end{bmatrix} = \begin{bmatrix}
x_D \\
y_D
\end{bmatrix} + t
\begin{bmatrix}
0 \\
1
\end{bmatrix}
$$
其中$t$是参数方程中的参数，$t>0$，以上向量表示的参数方程展开写就是
$$
\begin{cases}
x = x_D \\
y = y_D + t
\end{cases}
$$
对所有参数$t$大于0的情形，射线上的点满足以下坐标方程
$$
\begin{cases}
x = x_D \\
y \geq y_D
\end{cases}
$$
然后，我们写出直线AB的参数方程，设直线上的点是$(x,y)$
$$
\begin{bmatrix}
x \\
y
\end{bmatrix} = \begin{bmatrix}
x_A \\
y_A
\end{bmatrix} + t \begin{bmatrix}
x_B - x_A \\
y_B - y_A
\end{bmatrix}
$$
展开来写就是
$$
\begin{cases}
x = x_A + t (x_B - x_A) \\
y = y_A + t (y_B - y_A)
\end{cases}
$$
消去$t$，得到直角坐标方程
$$ \frac{x}{x_B-x_A} - \frac{y}{y_B-y_A} = \frac{x_A}{x_B-x_A} - \frac{y_A}{y_B-y_A} $$
$(x_A,y_A),(x_B,y_B)$是点A和点B的坐标，是已知的，而式中的$(x,y)$表示直线AB中的任何一点，不难理解，直线AB上的任意一点$(x,y)$是满足上面这个方程的.
同理，可得，直线BC的直角坐标方程
$$ \frac{x}{x_C-x_B} - \frac{y}{y_C-y_B} = \frac{x_B}{x_C-x_B} - \frac{y_B}{y_C-y_B} $$
和CA的直角坐标方程是
$$ \frac{x}{x_A-x_C} - \frac{y}{y_A-y_C} = \frac{x_C}{x_A-x_C} - \frac{y_C}{y_A-y_C} $$
我们实际上是想判断点D是不是在三角形ABC中，用射线法的思想，就是要判断从D出发的这条射线
$$
\begin{cases}
x = x_D \\
y \geq y_D
\end{cases}
$$
和直线AB，和直线BC，和直线CA这三条直线，也就是三角形ABC的三角形，总共相交多少次，如果相交的次数是奇数次，那么，点就是在三角形内部，如果，相交的次数是偶数次，那么，点就是在三角形外．具体理解，可以参考前面的那个动画演示.
「显然」，「不难」解得，如果有交点，如果从D朝着$\pi/2$方向出发的这条射线和AB，BC或者CD有交点，那么交点应该是：
$$ (x_D, y_1), (x_D, y_2), (x_D, y_3) $$
其中的$y_1$等于：
$$ y_1 = (\frac{x_D- x_A}{x_B-x_A} + \frac{y_A}{y_B-y_A}) \times (y_B - y_A) $$
$y_2$等于：
$$ y_2 = (\frac{x_D- x_B}{x_C-x_B} + \frac{y_B}{y_C-y_B}) \times (y_C - y_B) $$
$y_3$等于：
$$ y_3 = (\frac{x_D- x_C}{x_A-x_C} + \frac{y_C}{y_A-y_C}) \times (y_A - y_C) $$
我们接下来只需判断$y_1,y_2,y_3$中有多少个大于等于$y_D$就好啦！
randomPoint[] := RandomReal[{0, 1}, 2];
Clear[p1, p2, p3];
randomTriangle[] := Association[{
p1 -> randomPoint[],
p2 -> randomPoint[],
p3 -> randomPoint[]
}];
isPointXOutsideTriangleABC[pointX_, pointA_, pointB_, pointC_] := (
xa = pointA[[1]];
ya = pointA[[2]];
xb = pointB[[1]];
yb = pointB[[2]];
xc = pointC[[1]];
yc = pointC[[2]];
xd = pointX[[1]];
yd = pointX[[2]];
y1 = ((xd - xa)/(xb - xa) + (ya)/(yb - ya)) * (yb - ya);
y2 = ((xd - xb)/(xc - xb) + (yb)/(yc - yb)) * (yc - yb);
y3 = ((xd - xc)/(xa - xc) + (yc)/(ya - yc)) * (ya - yc);
counts = (Counts[{
y1 >= yd && xd >= Min[xa, xb] && xd <= Max[xa, xb],
y2 >= yd && xd >= Min[xb, xc] && xd <= Max[xb, xc],
y3 >= yd && xd >= Min[xc, xa] && xd <= Max[xc, xa]
}][True]) /. _Missing -> 0;
EvenQ[counts]
);
drawPointAndTriangle[pointX_, pointA_, pointB_, pointC_] := (
Graphics[
{
{
EdgeForm[{Thick, Gray}],
FaceForm[None],
Polygon[{pointA, pointB, pointC}]
},
{
PointSize[Large],
Point[pointX],
Gray
}
},
ImageSize -> {160, 160}
]
);
test[] := (
pX = randomPoint[];
rT = randomTriangle[];
pA = rT[p1];
pB = rT[p2];
pC = rT[p3];
{
drawPointAndTriangle[pX, pA, pB, pC],
isPointXOutsideTriangleABC[pX, pA, pB, pC]
}
)
samples = Table[test[], 500]
Counts[samples[[All, 2]]]
在Mathematica中检验我们的算法
上面我们写了一些代码，isPointXOutsideTriangleABC实际上包含了判定的过程代码，而random*函数呢生成随机点和随机三角形，drawPointAndTriangle把点和三角形画出来，好让我们亲自去看看这个点是不是在三角形中，以此查验程序给出的判定结果是不是正确的，上图是输出结果的一部．我自己查验了输出结果前面的200个，发现基本没什么问．值得注意的一点是，如果三角形的三个顶点是在$[0,1]$范围内随机取值，假设有一个随机的点$(x,y)$，其中$x \in [0,1]$，$y \in [0,1]$，那么$(x,y)$落入这个三角形的概率是很小的，Counts那行代码正是估算这个概率，怎么说呢，从实验结果估算，这个随机的点落入到这个随机的三角形中的概率，不到百分之8.
随机点落入随机三角形的概率估计
改天我们会再好好研究这个问题，这可以看做是这次的额外发现.
角度法
上面呢，我们已经展示了射线法基本思想，给出了射线法的数学推导过程，并且根据推导出的数学公式，写出了相应的这个用来判定一个点是不是在一个三角形内部的Mathematica程序代码，并且，还给出了程序输出结果的展．下面呢，我们会开始介绍另外一种，和射线法完全不一样的方法，角度法，顾名思义，是通过角度来判断一个点是否是在三角形内部．
我们来看一下上面这段演示，可以观察到呢，不管这个动点在三角形里边怎么动，它最大的那个角，都不会超过180度，这个动点越靠近中心，这些角度就越小，而最大的那个角度只有当这个点靠近三角形的边的时候，才会变．可以很容易地理解，当这个动点不小心移到三角形外边去的时候，它最大的那个角肯定会超过180度．于是乎，判定这个点是不是在三角形里边这个问题，就变成了，假设这个动点是P，判定，角APB，角BPC，角APC，这三个角中，有没有至少一个是大于180度的.
角度法的实现
我们得到的只是四个点的坐标，也就是那个点，还有那个三角形的三个点，分别记做$D, A, B, C．命题
$$D \in \triangle ABC$$
等价于
$$ \max(\angle ADB, \angle ADC, \angle CDB) \leq \pi $$
然后我们再把这四个点的坐标，分别记为$(x_D,y_D)$，$(x_A,y_A)$，$(x_B,y_B)$，$(x_C,y_C)．首先开始研究怎么样由坐标算出角度（或者至少是和角度相关的量来），我们假设有下面这样一个任意一个确定了坐标的三角形，三个顶点分别是P1，P2和P3，P1对边的斜率记做s1，P2对边的斜率记做s2，P3对边的斜率记做s3.
一个任意的三角形
三个顶点的坐标呢，分别记做$(x_{p1},y_{p1})$，$(x_{p2},y_{p2})$，和$(x_{p3},y_{p3})．假设我们要算P1这个顶点它对应的那个内角$\alpha_1 = \angle P_2 P_1 P_3$的余弦$\cos \alpha_1．过$P_2$做$P_1 P_3$的垂线，垂足为$P_2^{\prime}$，那么，根据余弦的定义，有
$$ \cos \alpha_1 = \frac{|P_1 P_2^{\prime}|}{|P_1 P_2|} $$
其中的$|P_1 P_2|$根据坐标，我们是可以很容易得出的，$|P_1 P_2^{\prime}|$则没那么显．不过，不妨先让我们假设，那条垂线$P_2 P_2^{\prime}$的方程是这样的斜截式，即，由截距和斜率确定的方程式
$$ y = \alpha_2^{\prime} + \beta_2^{\prime} x $$
其中的$(x,y)$表示这条垂线上的任意．下一步呢，首先，这个垂线的斜率可以很容易得到，两条互相垂直的直线的斜率互为负导数，因此，
$$ \beta_2^{\prime} = - \frac{1}{s_2}, \quad s_2 \neq 0 $$
如果在实际应用中$s_2$等于0，我们把这个三角形旋转一下就好．下一步要确定$\alpha_2^{\prime}$，现在已经有了斜率，并且知道这条直线经过的一点$P_2$，我们可以写出这条直线的参数方程
$$
\begin{bmatrix}
x \\
y
\end{bmatrix} = \begin{bmatrix}
x_{p2} \\
y_{p2}
\end{bmatrix} + t \begin{bmatrix}
1 \\
-s_2^{-1}
\end{bmatrix}
$$
要求截距，只需令方程中的$x$等于0，然后解出$y$，这里，在这个参数方程中，当$x$等于0时，$t$等于$-x_{p2}$，把$x=-x_{p2}$代入关于$y$的那部分式子中，得到
$$ y = y_{p2} + x_{p2} s_2^{-1} $$
这时的$y$就是截距，因为这是限定$x$取$0$才解出来的$y$，因此
$$ \alpha_2^{\prime} = y_{p2} + x_{p2} s_2^{-1} $$
现在，这条直线$P_2 P_2^{\prime}$的截斜式方程的两个参数$\alpha_2^{\prime}$以及$\beta_2^{\prime}$是基本确定．我们把直线$P_1 P_3$的方程记做
$$ y = \alpha_2 + \beta_2 x $$
注意，这和垂线的方程相比少了两个一撇上标，是两个方程，为了求交点，我们把这两个方程联立
$$
\begin{cases}
y = \alpha_2 + \beta_2 x \\
y = \alpha_2^{\prime} + \beta_2^{\prime}
\end{cases}
$$
改写这个方程组
$$
\begin{cases}
y - \beta_2 x = \alpha_2 \\
y - \beta_2^{\prime} x = \alpha_2^{\prime}
\end{cases}
$$
再把这个方程组写成矩阵形式
$$
\begin{bmatrix}
1, & -\beta_2 \\
1, & -\beta_2^{\prime}
\end{bmatrix} \begin{bmatrix}
y \\
x
\end{bmatrix}
= \begin{bmatrix}
\alpha_2 \\
\alpha_2^{\prime}
\end{bmatrix}
$$
注意，左边那个矩阵是2行2列的，比如说第一行，不是一减贝塔二，而是两个元素，一，和负的贝塔二，于是，交点亦可以用矩阵来表示，把这个交点记做$(x_{p2r}, y_{p2r})^{\mathsf{T}}$
$$
\begin{bmatrix}
y_{p2r} \\
x_{p2r}
\end{bmatrix} = \begin{bmatrix}
1, & -\beta_2 \\
1, & -\beta_2^{\prime}
\end{bmatrix}^{-1} \begin{bmatrix}
\alpha_2 \\
\alpha_2^{\prime}
\end{bmatrix}
$$
于是，
$$ \cos \alpha_1 = \frac{|P_1 P_2^{\prime}|}{|P_1 P_2|} = \frac{\sqrt{(x_{p2r}-x_1)^2+(y_{p2r}-y_1)^2}}{\sqrt{(x_2-x_1)^2+(y_2-y_1)^2}} $$
并且
$$ \sin \alpha_1 = \frac{|P_2 P_2^{\prime}|}{|P_1 P_2|} = \frac{\sqrt{(x_2-x_{2pr})^2+(y_2-y_{2pr})^2}}{\sqrt{(x_2-x_1)^2+(y_2-y_1)^2}} $$
现在回到原来的问题，请看图
点在三角形之外
我们看到，这时，即使点移出了三角形，仅仅根据坐标来算的话，得到的是内角$\alpha_1 + \alpha_2$而非我们真正所指的这里的外角$\alpha_3$，所以即使$\alpha_3$是大于$\pi$的，仅仅通过坐标计算出来的角度也看不．得变通一下，不能一直关注着坐标和角．我们还看到，点在三角形内部时，这三个角，加起来，显然是$2 \pi$，而点在外部时，$\alpha_3$的内角小于180度，而$\alpha_1 + \alpha_2$正是$\alpha_3$的内角，也小于180度，合着$\alpha_1$加上$\alpha_2$再加上$\alpha_3$的内角，得到的和角，小于$2 \pi$！这个结论有什么用呢？且看余弦函数的具有的对称性质
$$ \forall x \in \mathbb{R}, \cos (\pi - x) = \cos (\pi + x) $$
令$x = y - \pi$，就有
$$ \forall y \in \mathbb{R}, \cos y = \cos (2 \pi - y) $$
余弦函数性质
这意味着，我们可以验证是否有$\alpha_1 + \alpha_2 + \alpha_3$（这里的$\alpha_3$已经是开始指内角了）是否等于$2 \pi$，因为如果是的话，下列三个式子必定全成立，那就说明点是在内部：
$$
\begin{cases}
\cos \alpha_1 = \cos (\alpha_2 + \alpha_3) \\
\cos \alpha_2 = \cos (\alpha_1 + \alpha_3) \\
\cos \alpha_3 = \cos (\alpha_1 + \alpha_2)
\end{cases}
$$
如果不是的话，这三个式子必定至少有一个不成立，那就说明点是在外．思路豁然清晰起来！上列的式子等式的右端那些和角的余弦又怎么算呢？也很简单：
$$
\cos (\alpha_2 + \alpha_3) = \cos \alpha_2 \cos \alpha_3 - \sin \alpha_2 \sin \alpha_3 \\
\cos (\alpha_1 + \alpha_3) = \cos \alpha_1 \cos \alpha_3 - \sin \alpha_1 \sin \alpha_3 \\
\cos (\alpha_1 + \alpha_2) = \cos \alpha_1 \cos \alpha_2 - \sin \alpha_1 \sin \alpha_2
$$
在Mathematica中验证和角余弦公式
而我们已经知道怎样用坐标去表示$\cos \alpha$和$\sin \alpha$．下面，就可以开始编程实现了.
findIntersection[alpha1_, beta1_, alpha2_, beta2_] := (
sol = Inverse[
{
{1, -beta1},
{1, -beta2}
}
].{{alpha1}, {alpha2}};
Association[{y -> sol[[1]], x -> sol[[2]]}]
);
findVertical[pointA_, pointB_, pointC_] := (
slopeAB = (pointB[[2]] - pointA[[2]])/(pointB[[1]] - pointA[[1]]);
slopeAC = (pointC[[2]] - pointA[[2]])/(pointC[[1]] - pointA[[1]]);
slopeBC = (pointC[[2]] - pointB[[2]])/(pointC[[1]] - pointB[[1]]);
slopeABv = -1/slopeAB;
slopeACv = -1/slopeAC;
slopeBCv = -1/slopeBC;
intersectABv = pointC[[2]] + pointC[[1]]/slopeAB;
intersectACv = pointB[[2]] + pointB[[1]]/slopeAC;
intersectBCv = pointA[[2]] + pointA[[1]]/slopeBC;
Association[{
"slopeABv" -> slopeABv,
"slopeACv" -> slopeACv,
"slopeBCv" -> slopeBCv,
"intersectABv" -> intersectABv,
"intersectACv" -> intersectACv,
"intersectBCv" -> intersectBCv
}]
);
findVerticalFoot[pointA_, pointB_, pointC_] := (
slopeAB = (pointB[[2]] - pointA[[2]])/(pointB[[1]] - pointA[[1]]);
slopeAC = (pointC[[2]] - pointA[[2]])/(pointC[[1]] - pointA[[1]]);
slopeBC = (pointC[[2]] - pointB[[2]])/(pointC[[1]] - pointB[[1]]);
intersectAB = pointA[[2]] -
pointA[[1]]*(pointB[[2]] - pointA[[2]])/( pointB[[1]] - pointA[[1]]);
intersectAC = pointA[[2]] -
pointA[[1]]*(pointC[[2]] - pointA[[2]])/( pointC[[1]] - pointA[[1]]);
intersectBC = pointB[[2]] -
pointB[[1]]*(pointC[[2]] - pointB[[2]])/( pointC[[1]] - pointB[[1]]);
verticals = findVertical[pointA, pointB, pointC];
slopeABv = verticals["slopeABv"];
slopeACv = verticals["slopeACv"];
slopeBCv = verticals["slopeBCv"];
intersectABv = verticals["intersectABv"];
intersectACv = verticals["intersectACv"];
intersectBCv = verticals["intersectBCv"];
footAB = findIntersection[
intersectAB, slopeAB,
intersectABv, slopeABv
];
footAC = findIntersection[
intersectAC, slopeAC,
intersectACv, slopeACv
];
footBC = findIntersection[
intersectBC, slopeBC,
intersectBCv, slopeBCv
];
Association[{
"footAB" -> footAB,
"footAC" -> footAC,
"footBC" -> footBC
}]
);
twoPointDistance[pointX_, pointY_] := Sqrt[
(pointX[[1]] - pointY[[1]])^2 +
(pointX[[2]] - pointY[[2]])^2
];
下面我们试一下$A=(1,1)$，$B=(1.5,3)$，$C=(2,1.5)$这三个点，先在Mathematica中用我们写的函数算出$\triangle ABC$的垂足，如下：
在Mathematica中做出来的垂足坐标
然后我们在几何画板中做出这个三角形，用几何画板的命令做出垂足，显示的坐标如图
在几何画板中做出来的垂足坐标
可以看到，我们的Mathematica代码正确地计算出了$\triangle ABC$的垂．这时，有了垂足的坐标，理应可以开始计算余弦值了，但是事实上，仅仅从坐标，算出的实际上将会是余弦值的绝对值，因为算出来的距离总会是正的，所以，这里要定义「正角」和「负角」，因而也就有了正的余弦，和负的余弦：
正角
负角
如上图中的点IGJ组成的三角形$\triangle IGH$，其中的点$J$是垂足，假如说我们规定$\overrightarrow{GH}$的方向为正方向，那么，当$\angle IGH$大于0度小于90度时，$\overrightarrow{GJ}$的方向为正，线段$GJ$的长度可以「看做是」正的，而当$\angle IGH$大于90度小于180度时，$\overrightarrow{GJ}$的方向是与正方向$\overrightarrow{GH}$相反的，因此，可以把线段$GJ$的长度「看做是」负．正因此，你会看到对于$0 \leq x \leq \pi/2$，有$\cos x \geq 0$，对于$\pi/2 \leq x \leq \pi$，有$\cos x \leq 0$，余弦值的「正」和「负」就是这么来的.
所以，在我们的实际计算中，假如我们知道了$\triangle IGH$的三个点：点$I$，点$G$，和点$H$的坐标，以及垂足$J$的坐标，我们可以计算
$$ \cos \angle IGH = \frac{|GJ|}{|IG|} $$
并且，只需注意，当$|JH| > |GH|$时，使分子的$|GJ|$乘上一个负1以确保负值能传递到余弦中即．那么，接下来，我们就可以实现我们的根据三角形三个顶点坐标和三个垂足坐标计算余弦值的函数了：
computeCosineValueForTriangle[
pointA_,
pointB_,
pointC_,
verticalFootA_,
verticalFootB_,
verticalFootC_
] := (
signA = If[
Norm[
verticalFootB - pointC
] > Norm [
pointA - pointC
],
-1,
1
];
signB = If[
Norm[
verticalFootC - pointA
] > Norm [
pointB - pointA
],
-1,
1
];
signC = If[
Norm[
verticalFootA - pointB
] > Norm[
pointC - pointB
],
-1,
1
];
absCosineA = Divide[
Norm[verticalFootB - pointA],
Norm[pointA - pointB]
];
absCosineB = Divide[
Norm[verticalFootC - pointB],
Norm[pointB - pointC]
];
absCosineC = Divide[
Norm[verticalFootA - pointC],
Norm[pointA - pointC]
];
cosineA = absCosineA * signA;
cosineB = absCosineB * signB;
cosineC = absCosineC * signC;
Association[{
"cosineA" -> cosineA,
"cosineB" -> cosineB,
"cosineC" -> cosineC
}]
);
computeCosineValueForTriangle[
pointA_,
pointB_,
pointC_
] := (
foots = findVerticalFoot[pointA, pointB, pointC];
computeCosineValueForTriangle[
pointA,
pointB,
pointC,
{
foots["footBC"][x],
foots["footBC"][y]
},
{
foots["footAC"][x],
foots["footAC"][y]
},
{
foots["footAB"][x],
foots["footAB"][y]
}
]
);
尝试性地进行计算，得到的结果如图：
计算结果
另外，作为验证，在几何画板中把这三个角度量出来：
几何画板量出来的角度
接下来计算量出来的这三个角度，看得到的那个输出结果是否正确：
比较两个结果
可以看到，我们算出来的余弦值，和直接量出来的余弦值，两者之间的差别是比较小的，可以认为我们是算对．至于这个差别是怎么来的呢？在几何画板中量这个角度的时候，只保留到小数点后一位，我估计，这两组数据的误差就是这么来的，我们是直接从坐标开始计算的，而在几何画板中，从坐标到量出来的角度有一个误差，造成了信息的损失，所以，可以说，我们的更精确一点.
接下来计算正弦值，正弦值就好说了，在$0 \leq x \leq \pi$都有$\cos x \geq 0$，所以，很简单，直接套用下面这个公式就可以了
$$ \sin x = \sqrt{1 - (\cos x)^2 } $$
computeSineValueForTriangle[pointA_, pointB_, pointC_] := (
cosines = computeCosineValueForTriangle[pointA, pointB, pointC];
Association[{
"sineA" -> Sqrt[1 - (cosines["cosineA"])^2],
"sineB" -> Sqrt[1 - (cosines["cosineB"])^2],
"sineC" -> Sqrt[1 - (cosines["cosineC"])^2]
}]
);
然后再来验证看：
比较和验证
误差没有太大，余弦值是对的，因此正弦值应该也是对的.
再下边，我们可以回到最开始的时候问的问题，给定一个点，叫做$D$，给定三个坐标$A$，$B$，和$C$，判断点$D$是否在$\triangle ABC$．而在开始计算余弦值和正弦值之前，我们已经有结论，$D \in \triangle ABC$ 当且仅当以下式子同时成立
$$
\begin{cases}
\cos \angle ADB = \cos(\angle BDC + \angle ADC) \\
\cos \angle BDC = \cos(\angle ADB + \angle ADC) \\
\cos \angle ADC = \cos(\angle ADB + \angle BDC)
\end{cases}
$$
其中，这四个点的位置如图：
四个点的位置
这样就可以开始实现了：
cosineAdd[cosineAlpha_, cosBeta_] := (
cosineAlpha * cosBeta -
Sqrt[1 - cosineAlpha^2] * Sqrt[1 - cosBeta^2]
);
tellIfPointDInsideTriangleABC[
cosineADB_,
cosineBDC_,
cosineADC_
] := (
epsilon = 10^(-1);
FullSimplify[
Abs[cosineADB - cosineAdd[cosineBDC, cosineADC]] < epsilon &&
Abs[cosineBDC - cosineAdd[cosineADB, cosineADC]] < epsilon &&
Abs[cosineADC - cosineAdd[cosineADB, cosineBDC]] < epsilon
]
);
tellIfPointDInsideTriangleABC[
pointD_,
pointA_,
pointB_,
pointC_
] := (
cosineADB = computeCosineValueForTriangle[
pointD, pointA, pointB
]["cosineA"];
cosineBDC = computeCosineValueForTriangle[
pointD, pointB, pointC
]["cosineA"];
cosineADC = computeCosineValueForTriangle[
pointD, pointA, pointC
]["cosineA"];
tellIfPointDInsideTriangleABC[
cosineADB,
cosineBDC,
cosineADC
]
);
接下来我们进行一下验证
这些点的位置
验证
可以看到，针对这个样例，咱们的程序给出了正确的结．角度法我这里给出的并不是一个好的实现，仅仅是为了验证它是理论上可行的，没有好的实现步骤，实现起来真的好麻烦，虽然原理都很简单都很容易理解.
面积法
要用面积法来判断一个点是不是在多边形内部其实也非常简单，假设知道了这么一个点$D$，点$A$，点$B$，以及点$C$确定的三角形$\triangle ABC$，那么，如果点$D$是在三角形内部，那么，点$D$和三个顶点的连线其实可以把整个三角形$\triangle ABC$分成三部分
$$ \triangle ABC = \triangle ABD + \triangle BDC + \triangle ADC $$
当点在三角形内部时，这个公式成立，如下图所演示的
并且总三角形的面积，实际上也等于三个小三角形的面积之和，只要点$D$在三角形$\triangle ABC$中
$$ S_{\triangle ABC} = S_{\triangle ABD} + S_{\triangle BDC} + S_{\triangle ADC} $$
那假如说点$D$在三角形$\triangle ABC$之外呢，这个面积等式还成立吗？
当点D移出三角形外
当我们试着把点$D$移出三角形$\triangle ABC$外之和，发现那个面积等式很明显就不成立．于是，很简单，要判断一个点$D$是不是在三角形$\triangle ABC$中，只需判断关系式
$$ S_{\triangle ABC} = S_{\triangle ABD} + S_{\triangle BDC} + S_{\triangle ADC} $$
是否成立就可以了.
面积法的实现
我们要实现这样一个判定函数，它接受四个坐标，当第一个坐标是位于后三个坐标所确定的三角形内部的时候，返回真值，当第一个坐标是位于后三个坐标所确定的三角形外部的时候，返回假．现在我们只知道点的坐标，而要用面积法判断的话，要知道三角形的面积，面积可以，根据海伦-秦九韶公式，用三角形的边长求出，而三角形的边长，则可以用勾股定理，由坐标求．
海伦-秦九韶公式是这样子的，设面积为$S_{\triangle}$，$a,b,c$分别是边长，那么有
$$ S_{\triangle} = \sqrt{s(s-a)(s-b)(s-c)} $$
式中的$s$，它等于
$$ s = \frac{a+b+c}{2} $$
注意区分这里的$S$的大小写.
下面首先实现计算面积的程序
computeTriangleArea[a_, b_, c_] := (
s = (a+b+c)/2;
Sqrt[s*(s-a)*(s-b)*(s-c)]
)
computeTriangleArea[
{xa_, ya_},
{xb_, yb_},
{xc_, yc_}
] := (
computeTriangleArea[
Norm[{xa-xb, ya-yb}],
Norm[{xa-xc, ya-yc}],
Norm[{xb-xc, yb-yc}]
]
);
非常简短，我们就不验证．接下来实现判断的部分
tellIfPointDInsideTriangleABC[
pointD_,
pointA_,
pointB_,
pointC_
] := (
epsilon = 10^(-1);
Abs[
computeTriangleArea[
pointA,
pointB,
pointC
] - computeTriangleArea[
pointA,
pointB,
pointD
] - computeTriangleArea[
pointB,
pointD,
pointC
] - computeTriangleArea[
pointA,
pointD,
pointC
]
] < epsilon
)
验证一下，要用来验证的点的位置是这样的：
点的位置
验证结果是这样的：
验证结果
面积法容易实现，容易理解，并且不容易出错，在二维空间中，面积法是很好用的.
代数法
在高中的时候有人可能已经学过了一些初步的线性规划方面的知识，一般叫做最优化问题，具体是，例如，有这样的问题，生产一件a商品能获得利润$p_1$元，生产一件b商品能获得利润$p_2$元，生产一件a商品消耗原材料甲$c_{11}$千克，原材料乙$c_{12}$千克，生产一件b商品消耗原材料甲$c_{21}$千克，原材料乙$c_{22}$千克，现有原材料甲$q_1$千克，原材料乙$q_2$千克，问你，要生产a商品多少件，b商品多少件，能获得最大利润？
这样的问题呢，是一个典型的线性规划问题，可以设生产$x$件商品a，$y$件商品b，那么，该问题可用线性规划模型来建模
$$
\text{maximize} \quad p_1 x + p_2 y \\
$$
$$
\text{subjected to} \quad \begin{cases}
c_{11} x + c_{21} y \leq q_1 \\
c_{12} x + c_{22} y \leq q_2 \\
x, y \geq 0
\end{cases}
$$
说了这么多，它跟判断一个点在不在三角形里面这个问题又有什么关系呢？先回想一下当时我们做这类题目的一种解法，我们具体是，把约束条件里的不等式改写成等式然后统统画在坐标系上，在被包围的那块区域找最优解，例如，画出来可能会是这样的：
可行域
约束条件其实是四条不等式，我们把不等号全部换做等号，就得到四条直线，例如，把$x \geq 0$这个不等式改写成$x = 0$这个等式，它就表示一条过原点垂直$x$轴的直线，而显然，$x \geq 0$描述的是一个区域，具体是，这条直线，这条过原点垂直$x$轴的直线右边那部分的区域，同理，$y \geq 0$描述的是过原点垂直$y$轴的那条直线上边的区域，而$c_{11}x + c_{21}y = q_1$和$c_{12}x+c_{22}y = q_2$应该代表的是两条斜线，不等式$c_{11}x + c_{21}y \leq q_1$和$c_{12}x+c_{22}y \leq q_2$表示的应该就是同时位于这两条斜线下方的区域，那么，这四个区域的公共交集，观察上面那幅图，你应该已经知道在哪了.
看到了吗？也是一个多边．即，一个多边形其实是可以用若干条不等式来表述的，因为一条不等式描述的是空间中的一个区域，而，若干条不等式，描述的仍然是空间中的区域只不过该区域同时是多个区域的交．由此我们立刻受到启发，是否一个三角形，也能用三条不等式来表示呢？如果是的话，那三角形内部的点，实际上就应该是同时使那三条不等式同时满足的点，而如果是三角形外部的点的话，那肯定会至少违背其中的一条不等式，这样不就可以解决我们要解决的问题了吗？
其实，平面上的三角形，是可以用三条不等式来表示的.
代数法的实现
如果要应用代数法的思想来解决我们要解决的问题，首先我们就要把三角形的三个坐标，转化为平面直角坐标系上描述这个三角形的三个不等式，具体地，不妨设三角形$\triangle ABC$的三个顶点分别为点$A$，点$B$，和$C$，三个顶点的坐标分别为$(x_A, y_A)$，$(x_B, y_B)$，和$(x_C, y_C)$，回忆我们关于「射线法」的讨论，可以立刻写出三角形三条边的直线方程：
$$
\frac{x}{x_B-x_A} - \frac{y}{y_B-y_A} = \frac{x_A}{x_B-x_A} - \frac{y_A}{y_B-y_A} \\
\frac{x}{x_C-x_B} - \frac{y}{y_C-y_B} = \frac{x_B}{x_C-x_B} - \frac{y_B}{y_C-y_B} \\
\frac{x}{x_A-x_C} - \frac{y}{y_A-y_C} = \frac{x_C}{x_A-x_C} - \frac{y_C}{y_A-y_C}
$$
现在问题是，要把三角形三条边的直线方程改写成不等式，那就要知道不等式的不等号的开口是朝哪边的，那是朝哪边的呢？我们可以利用三角形的重心所具有的性质，一个三角形的重心，不管这个三角形是什么形状的，这个重心总是会在三角形的内部而不会是在外部，我们现在只需利用这个性质即可，过后再去证明．重心既然总是在里边的，那么，如果这三个不等式的不等号设置得正确的话，重心就会同时满足这三条不等式，我们逐一尝试所有可能的不等号的配置组合直到找到那个正确的组合为止即可.
当然了，首先，要把这个重心的位置给它找出来，其实不是很复杂，有公式，重心的坐标是
$$ (\frac{1}{3}(x_A + x_B + x_C), \frac{1}{3}(y_A + y_B + y_C)) $$
顶点坐标我们是知道的，带进去算，就可以得到重心的坐标．下面开始判断
tellIfPointDInsideTriangleABC[
{xd_, yd_},
{xa_, ya_},
{xb_, yb_},
{xc_, yc_}
] := (
xm = (xa+xb+xc)/3;
ym = (ya+yb+yc)/3;
test1 = FullSimplify[
xm/(xb-xa)-ym/(yb-ya) > xa/(xb-xa)-ya/(yb-ya) &&
xm/(xc-xb)-ym/(yc-yb) < xb/(xc-xb)-yb/(yc-yb) &&
xm/(xa-xc)-ym/(ya-yc) < xc/(xa-xc)-yc/(ya-yc) &&
xd/(xb-xa)-yd/(yb-ya) > xa/(xb-xa)-ya/(yb-ya) &&
xd/(xc-xb)-yd/(yc-yb) < xb/(xc-xb)-yb/(yc-yb) &&
xd/(xa-xc)-yd/(ya-yc) < xc/(xa-xc)-yc/(ya-yc)
];
test2 = FullSimplify[
xm/(xb-xa)-ym/(yb-ya) < xa/(xb-xa)-ya/(yb-ya) &&
xm/(xc-xb)-ym/(yc-yb) > xb/(xc-xb)-yb/(yc-yb) &&
xm/(xa-xc)-ym/(ya-yc) < xc/(xa-xc)-yc/(ya-yc) &&
xd/(xb-xa)-yd/(yb-ya) < xa/(xb-xa)-ya/(yb-ya) &&
xd/(xc-xb)-yd/(yc-yb) > xb/(xc-xb)-yb/(yc-yb) &&
xd/(xa-xc)-yd/(ya-yc) < xc/(xa-xc)-yc/(ya-yc)
];
test3 = FullSimplify[
xm/(xb-xa)-ym/(yb-ya) < xa/(xb-xa)-ya/(yb-ya) &&
xm/(xc-xb)-ym/(yc-yb) < xb/(xc-xb)-yb/(yc-yb) &&
xm/(xa-xc)-ym/(ya-yc) > xc/(xa-xc)-yc/(ya-yc) &&
xd/(xb-xa)-yd/(yb-ya) < xa/(xb-xa)-ya/(yb-ya) &&
xd/(xc-xb)-yd/(yc-yb) < xb/(xc-xb)-yb/(yc-yb) &&
xd/(xa-xc)-yd/(ya-yc) > xc/(xa-xc)-yc/(ya-yc)
];
test4 = FullSimplify[
xm/(xb-xa)-ym/(yb-ya) < xa/(xb-xa)-ya/(yb-ya) &&
xm/(xc-xb)-ym/(yc-yb) > xb/(xc-xb)-yb/(yc-yb) &&
xm/(xa-xc)-ym/(ya-yc) > xc/(xa-xc)-yc/(ya-yc) &&
xd/(xb-xa)-yd/(yb-ya) < xa/(xb-xa)-ya/(yb-ya) &&
xd/(xc-xb)-yd/(yc-yb) > xb/(xc-xb)-yb/(yc-yb) &&
xd/(xa-xc)-yd/(ya-yc) > xc/(xa-xc)-yc/(ya-yc)
];
test5 = FullSimplify[
xm/(xb-xa)-ym/(yb-ya) > xa/(xb-xa)-ya/(yb-ya) &&
xm/(xc-xb)-ym/(yc-yb) < xb/(xc-xb)-yb/(yc-yb) &&
xm/(xa-xc)-ym/(ya-yc) > xc/(xa-xc)-yc/(ya-yc) &&
xd/(xb-xa)-yd/(yb-ya) > xa/(xb-xa)-ya/(yb-ya) &&
xd/(xc-xb)-yd/(yc-yb) < xb/(xc-xb)-yb/(yc-yb) &&
xd/(xa-xc)-yd/(ya-yc) > xc/(xa-xc)-yc/(ya-yc)
];
test6 = FullSimplify[
xm/(xb-xa)-ym/(yb-ya) > xa/(xb-xa)-ya/(yb-ya) &&
xm/(xc-xb)-ym/(yc-yb) > xb/(xc-xb)-yb/(yc-yb) &&
xm/(xa-xc)-ym/(ya-yc) < xc/(xa-xc)-yc/(ya-yc) &&
xd/(xb-xa)-yd/(yb-ya) > xa/(xb-xa)-ya/(yb-ya) &&
xd/(xc-xb)-yd/(yc-yb) > xb/(xc-xb)-yb/(yc-yb) &&
xd/(xa-xc)-yd/(ya-yc) < xc/(xa-xc)-yc/(ya-yc)
];
(test1) ||
(test2) ||
(test3) ||
(test4) ||
(test5) ||
(test6)
);
下面开始验证，这是验证用的图形
代数法用来验证的图形
这是验证结果.
代数法的验证
代数法只用一个函数就能够搞得定，是符合「高内聚低耦合」编程理念的典范，看起来，和面积法一样容．并且，实际上你会发现，代数法可以轻松地推广到更高维的空间，也就是说，不仅是三角形，就算是一般的n边形，甚至，3维空间中的n面体，都不在话．我们这里就不一一实现了，这个，交给感兴趣的读．
按照难度来比较的话，我觉得最容易的是面积法，然后是代数法，然后是射线法，最麻烦的是角度法，最好不要去用角度法来判断点是不是在三角形里面，因为复杂，而且容易出错.
总结
我们这篇文章一开始提出了这样一个问题：给定一个点，给定一个多边形，怎样根据点和多边形顶点的坐标来判断这个点是不是在多边形内部，后来，我们为了简化问题，把讨论的范围限制在2维空间和多边形的边数为3（三角形）的情．我们先后介绍了射线法、角度法、面积法、代数法总共这四种方法，其中的每一种方法都能判断一个点是否是在三角形内部，射线法的思路是：点在移动时，其穿入多边形边界的次数等于其穿出多边形次数的边界，因而如果点在外部，穿入边界次数等于穿出边界次数，射线与多边形的边的总交点数是偶数，角度发的思路是：点在三角形内部时，点到三角形顶点的连线夹成的3个角的和总是360度，而点在外部时，这三个角的和则小于360度，和角余弦公式可以帮助验证这一点，面积法的思路是：点在三角形内部时，点到三个顶点的连线和原来三角形的三条边总共把这个三角形切成三部分，而这三部分的面积和，等于原来的三角形，点在外部时，这三部分的面积和，大于原来的三角形，而有了坐标，我们可以用海伦公式算出面积，通过比较面积，判断点是否在三角形内部，代数法的思路是：把三角形区域看成是3个不等式描述的范围的交集，如果点在这个区域内，那么点必定满足那三个不等式，把点代入不等式去算就知道．
这样一个看似繁琐无聊的问题，实则有很大的用．在游戏的逻辑判断模块中，需要知道一个角色是不是已经到达了地图边界，而整张地图实际上就是一个平面上的多边形，需要在角色到达边界时限制其再向外移动，在物理引擎的实现，或者仿真模拟的过程中，可能会要模拟到粒子的碰撞，一个粒子可以看做是多面体，要判断两个粒子是否已经撞在了一起，需要判断其中一个多面体的顶点有没有在另外一个多面体内部，这是3维空间中的覆盖判定问题.
我们这里给出的所谓「实现」，与其说是实现，不如说是一种「验证」，即，验证我们的推理是不是正确的，Mathematica具有简单优雅的语法，和强大的数学功能，因此，我经常会选用Mathematica来验证自己的猜想，当然了，推理的正确性不能全靠验证，我们事实上也是真正地进行过分析和讨论来着的，公式和演算也尽可能地做得严谨，所以，如果您看了觉得这些个公式的推导过程是没问题的，大概率就是没问题的，如果有问题，也可以跟我．Mathematica的性能，和真正的生产环境是差得远的，计算过程是类似的，你可以尝试着在Python中实现这些过程，我不是太喜欢C++，所以如果不是被逼迫，应该不会写C++的，你可以在任何语言中尝试实现这些计算和判定的过程.
交流与探讨
很有可能，判断一个点是否在多边形内部还有好多种方法，但是我现在知道的，暂时也就这几种，如果你知道别的更简单，更高效的方法，可以告诉我.
如果，你实现了点与矩形，更一般的四边形，甚至一般的平面n边形，甚至空间或者更高维空间的n面体的覆盖判定问题的判定函数，我很期待你和我分享思路、数学推导过程甚至程序代码（不一定要是Mathematica代码），当然，有详细的数学推导过程最好.
我的邮箱是
[email protected]
可以用Wayne这个署名来称呼我，如有问题探讨和交流、好主意分享、答疑、质询、勘误、友链交换等目的都可以通过这个邮箱与我联．
感谢你阅读到这，如果你觉得我写得对你有帮助，欢迎订阅这个网站的RSS，或者如果你也经营着自己的博客网站，可以和我申请友链.